search

SciML / ModelingToolkit.jl

An acausal modeling framework for automatically parallelized scientific machine learning (SciML) in Julia. A computer algebra system for integrated symbolics for physics-informed machine learning and automated transformations of differential equations
https://mtk.sciml.ai/dev/
Other
1.43k stars 209 forks source link

Compose Vector and Scalar Equations together #2986

Open Snowd1n opened 2 months ago

Snowd1n commented 2 months ago

When I try to connect an equation of scalars with an equation of vectors, something on the compose step fails even though what is being connected is the same type (both scalar etc.). I am not sure why this is. An MWE is below

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

function FOL(;name)
    @parameters begin
        τ[1:5]  # parameters
    end
    @variables begin
        x(t)[1:5] = 0.0
        a(t) # dependent variables
    end

    eqs = D(x) ~ (a .- x) ./ τ

    ODESystem(eqs,t;name)
end

function FOL_scal(;name)
    @parameters begin
        b = 0.5
        c = 4.0 # parameters
    end
    @variables begin
        a(t) = 1.0 # dependent variables
    end

    eqs = D(a) ~ (a - b)/c

    ODESystem(eqs,t;name)
end

@named fol = FOL()
@named scal = FOL_scal()
connected = compose(ODESystem([scal.a ~ fol.a],t,name=:connected,defaults=Dict(),discrete_events=[]),fol,scal)

Error & Stacktrace ⚠️

ERROR: MethodError: no method matching length(::SymbolicUtils.BasicSymbolic{Any})

Closest candidates are:
  length(::Tables.DictRowTable)
   @ Tables C:\Users\u5522838\.julia\packages\Tables\8p03y\src\dicts.jl:118
  length(::CSTParser.EXPR)
   @ CSTParser C:\Users\u5522838\.julia\packages\CSTParser\0hXvH\src\spec.jl:278
  length(::Cmd)
   @ Base process.jl:678
  ...

Stacktrace:
  [1] _similar_shape(itr::SymbolicUtils.BasicSymbolic{Any}, ::Base.HasLength)
    @ Base .\array.jl:710
  [2] _collect(cont::UnitRange{Int64}, itr::SymbolicUtils.BasicSymbolic{Any}, ::Base.HasEltype, isz::Base.HasLength)
    @ Base .\array.jl:765
  [3] collect(itr::SymbolicUtils.BasicSymbolic{Any})
    @ Base .\array.jl:759
  [4] broadcastable(x::SymbolicUtils.BasicSymbolic{Any})
    @ Base.Broadcast .\broadcast.jl:743
  [5] broadcasted
    @ .\broadcast.jl:1344 [inlined]
  [6] broadcast(::typeof(-), ::SymbolicUtils.BasicSymbolic{Any}, ::SymbolicUtils.BasicSymbolic{Vector{Real}})
    @ Base.Broadcast .\broadcast.jl:841
  [7] maketerm(::Type{Symbolics.ArrayOp{AbstractVector{Real}}}, f::Function, args::Vector{Any}, _symtype::Type, m::Nothing)
    @ Symbolics C:\Users\u5522838\.julia\packages\Symbolics\qKoME\src\arrays.jl:66
  [8] namespace_expr(O::Symbolics.ArrayOp{AbstractVector{Real}}, sys::ODESystem, n::Symbol; ivs::Vector{SymbolicUtils.BasicSymbolic{Real}})
    @ ModelingToolkit C:\Users\u5522838\.julia\packages\ModelingToolkit\GJiqn\src\systems\abstractsystem.jl:1269
  [9] namespace_expr
    @ C:\Users\u5522838\.julia\packages\ModelingToolkit\GJiqn\src\systems\abstractsystem.jl:1250 [inlined]
 [10] (::ModelingToolkit.var"#298#300"{Vector{SymbolicUtils.BasicSymbolic{Real}}, ODESystem, Symbol})(a::Symbolics.ArrayOp{AbstractVector{Real}})
    @ ModelingToolkit C:\Users\u5522838\.julia\packages\ModelingToolkit\GJiqn\src\systems\abstractsystem.jl:1258
 [11] iterate
    @ .\generator.jl:47 [inlined]
 [12] collect_to!(dest::Vector{typeof(/)}, itr::Base.Generator{Vector{Any}, ModelingToolkit.var"#298#300"{Vector{SymbolicUtils.BasicSymbolic{Real}}, ODESystem, Symbol}}, offs::Int64, st::Int64)
    @ Base .\array.jl:892
 [13] collect_to_with_first!(dest::Vector{typeof(/)}, v1::Function, itr::Base.Generator{Vector{Any}, ModelingToolkit.var"#298#300"{Vector{SymbolicUtils.BasicSymbolic{Real}}, ODESystem, Symbol}}, st::Int64)
    @ Base .\array.jl:870
 [14] _collect(c::Vector{Any}, itr::Base.Generator{Vector{Any}, ModelingToolkit.var"#298#300"{Vector{SymbolicUtils.BasicSymbolic{Real}}, ODESystem, Symbol}}, ::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base .\array.jl:864
 [15] collect_similar(cont::Vector{Any}, itr::Base.Generator{Vector{Any}, ModelingToolkit.var"#298#300"{Vector{SymbolicUtils.BasicSymbolic{Real}}, ODESystem, Symbol}})
    @ Base .\array.jl:763
 [16] map(f::Function, A::Vector{Any})
    @ Base .\abstractarray.jl:3285
 [17] namespace_expr(O::Symbolics.ArrayOp{AbstractVector{Real}}, sys::ODESystem, n::Symbol; ivs::Vector{SymbolicUtils.BasicSymbolic{Real}})
    @ ModelingToolkit C:\Users\u5522838\.julia\packages\ModelingToolkit\GJiqn\src\systems\abstractsystem.jl:1258
 [18] namespace_equation(eq::Equation, sys::ODESystem, n::Symbol; ivs::Vector{SymbolicUtils.BasicSymbolic{Real}})
    @ ModelingToolkit C:\Users\u5522838\.julia\packages\ModelingToolkit\GJiqn\src\systems\abstractsystem.jl:1240
 [19] namespace_equation (repeats 2 times)
    @ C:\Users\u5522838\.julia\packages\ModelingToolkit\GJiqn\src\systems\abstractsystem.jl:1235 [inlined]
 [20] #292
    @ C:\Users\u5522838\.julia\packages\ModelingToolkit\GJiqn\src\systems\abstractsystem.jl:1225 [inlined]
 [21] iterate
    @ .\generator.jl:47 [inlined]
 [22] _collect(c::Vector{Equation}, itr::Base.Generator{Vector{Equation}, ModelingToolkit.var"#292#293"{ODESystem, Vector{SymbolicUtils.BasicSymbolic{Real}}}}, ::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base .\array.jl:854
 [23] collect_similar(cont::Vector{Equation}, itr::Base.Generator{Vector{Equation}, ModelingToolkit.var"#292#293"{ODESystem, Vector{SymbolicUtils.BasicSymbolic{Real}}}})
    @ Base .\array.jl:763
 [24] map(f::Function, A::Vector{Equation})
    @ Base .\abstractarray.jl:3285
 [25] namespace_equations(sys::ODESystem, ivs::Vector{SymbolicUtils.BasicSymbolic{Real}})
    @ ModelingToolkit C:\Users\u5522838\.julia\packages\ModelingToolkit\GJiqn\src\systems\abstractsystem.jl:1225
 [26] namespace_equations(sys::ODESystem)
    @ ModelingToolkit C:\Users\u5522838\.julia\packages\ModelingToolkit\GJiqn\src\systems\abstractsystem.jl:1223
 [27] _broadcast_getindex_evalf
    @ .\broadcast.jl:709 [inlined]
 [28] _broadcast_getindex
    @ .\broadcast.jl:682 [inlined]
 [29] getindex(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(ModelingToolkit.namespace_equations), Tuple{Base.Broadcast.Extruded{Vector{ODESystem}, Tuple{Bool}, Tuple{Int64}}}}, I::Int64)
    @ Base.Broadcast .\broadcast.jl:636
 [30] copy(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(ModelingToolkit.namespace_equations), Tuple{Vector{ODESystem}}})
    @ Base.Broadcast .\broadcast.jl:942
 [31] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(ModelingToolkit.namespace_equations), Tuple{Vector{ODESystem}}})
    @ Base.Broadcast .\broadcast.jl:903
 [32] equations(sys::ODESystem)
    @ ModelingToolkit C:\Users\u5522838\.julia\packages\ModelingToolkit\GJiqn\src\systems\abstractsystem.jl:1464
 [33] show(io::IOContext{Base.TTY}, mime::MIME{Symbol("text/plain")}, sys::ODESystem)
    @ ModelingToolkit C:\Users\u5522838\.julia\packages\ModelingToolkit\GJiqn\src\systems\abstractsystem.jl:1770
 [34] (::REPL.var"#55#56"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL C:\Users\u5522838\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:273
 [35] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
    @ REPL C:\Users\u5522838\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:569
 [36] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL C:\Users\u5522838\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:259
 [37] display(d::REPL.REPLDisplay, x::Any)
    @ REPL C:\Users\u5522838\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:278
 [38] display(x::Any)
    @ Base.Multimedia .\multimedia.jl:340
 [39] #invokelatest#2
    @ .\essentials.jl:892 [inlined]
 [40] invokelatest
    @ .\essentials.jl:889 [inlined]
 [41] (::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\u5522838\.vscode\extensions\julialang.language-julia-1.83.2\scripts\packages\VSCodeServer\src\eval.jl:237
 [42] withpath(f::VSCodeServer.var"#69#74"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
    @ VSCodeServer c:\Users\u5522838\.vscode\extensions\julialang.language-julia-1.83.2\scripts\packages\VSCodeServer\src\repl.jl:276
 [43] (::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\u5522838\.vscode\extensions\julialang.language-julia-1.83.2\scripts\packages\VSCodeServer\src\eval.jl:179
 [44] hideprompt(f::VSCodeServer.var"#68#73"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
    @ VSCodeServer c:\Users\u5522838\.vscode\extensions\julialang.language-julia-1.83.2\scripts\packages\VSCodeServer\src\repl.jl:38
 [45] (::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\u5522838\.vscode\extensions\julialang.language-julia-1.83.2\scripts\packages\VSCodeServer\src\eval.jl:150
 [46] with_logstate(f::Function, logstate::Any)
    @ Base.CoreLogging .\logging.jl:515
 [47] with_logger
    @ .\logging.jl:627 [inlined]
 [48] (::VSCodeServer.var"#66#71"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer c:\Users\u5522838\.vscode\extensions\julialang.language-julia-1.83.2\scripts\packages\VSCodeServer\src\eval.jl:263
 [49] #invokelatest#2
    @ .\essentials.jl:892 [inlined]
 [50] invokelatest(::Any)
    @ Base .\essentials.jl:889

Version 9.32.0 of ModelingToolkit and Version 1.10.4 of Julia

Additional context

Add any other context about the problem here.

Snowd1n commented 2 months ago

Contradict on Discourse suggested using a semicolon at the end of the line as it appears that the REPL can't display the connected equation. However this then fails at tstructural_simplify(connected) with the same top level error. THe URL to that discussion is here: https://discourse.julialang.org/t/compose-vector-and-scalar-equations-together/118474

Snowd1n commented 2 months ago

Reply from Contradict

OK, thats a different stack trace and a different problem. I am no longer sure if either of these is a bug, but you can fix both by broadcasting the equation in FOL!

using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D using DifferentialEquations

function FOL(; name) @parameters begin τ[1:5] # parameters end @variables begin x(t)[1:5] = 0.0 a(t) # dependent variables end

eqs = D(x) .~ (a .- x) ./ τ # Changing this to .~

ODESystem(eqs, t; name)

end

function FOL_scal(; name) @parameters begin b = 0.5 c = 4.0 # parameters end @variables begin a(t) = 1.0 # dependent variables end

eqs = D(a) ~ (a - b) / c

ODESystem(eqs, t; name)

end

function test() @named fol = FOL() @named scal = FOL_scal() connected = compose( ODESystem( [scal.a ~ fol.a], t, name = :connected, defaults = Dict(), discrete_events = [], ), fol, scal, )

simp = complete(structural_simplify(connected)) # complete makes variable and parameter access easier.

prob = ODEProblem(simp, [], (0.0, 10.0), [simp.fol.τ=>ones(5)]) # I added this prameter value too

sol = solve(prob, Tsit5(); abstol = 1e-3, reltol = 1e-3)

end The reason I am no longer sure this is a bug is because I have seen other cases where removing all the broadcast operations and leaving equations as all vectors works too. When you get stack traces involving the broadcast machinery, it is worth trying both to see if one works.

Snowd1n commented 2 months ago

That now works I agree. But in my real case, this doesn’t solve the problem and the broadcasting creates the Array has no term lhs error which I know from other threads is fixed via “eqs = Symbolics.scalarize(reduce(vcat, Symbolics.scalarize.(eqs)))”. However, this then leads back to the same structural_simplify error as before. I don’t really know where to proceed from here as I basically have now tried both variants and loop to the same point. I think it may be linked to scalarize