Compose Vector and Scalar Equations together

Specific Domains Modelling & Simulations
, , ,
1

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.

An issue has been raised on the ModelingToolkit git : Compose Vector and Scalar Equations together · Issue #2986 · SciML/ModelingToolkit.jl · GitHub.
But I wanted to post here for anyone else to see or help as I think this may be me doing something wrong. I shall add any comments on either to the respective forum.

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: 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
2

This does seem like a bug, but I think it is in the show method. Try adding a semicolon to that last line, connected may actually be ok, you just can’t print it:

connected = compose(ODESystem([scal.a ~ fol.a],t,name=:connected,defaults=Dict(),discrete_events=[]),fol,scal);
3

Unfortunately, I don’t think that’s the case as the structural simplify command then fails on this connected system and I can’t get past that with another semicolon. Ngl, I don’t really know why the semicolon makes that line work but it still sticks at the next.

simp = structural_simplify(connected); #Fails here 

prob=ODEProblem(simp,[],(0.0,10.0),[]);

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

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.

5

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.

6

Yeah, that sounds like the same battle. Can you post a more complete example? I usually avoid manually scalarizing unless absolutely necessary, most of the time the problem can be solved by broadcasting everything or broadcasting nothing.

7

Hello, thank you for all your help. Here is a more complex example. Basically I need the function athem_factory to be able to be simplified but I can’t make that happen.

using ModelingToolkit,Symbolics,Dierckx,Integrals
using ModelingToolkit: t_nounits as t, D_nounits as D

function athem_factory(DOS::Spline1D,laser::Num,egl::Int;name)
    @parameters Tel

    @named dfneq = athem_template(egl)

    connections = [dfneq.Source ~ athem_excitation(dfneq.fneq,Tel,laser,DOS,egl)]

    compose(ODESystem(connections,t;name),dfneq)
end

function athem_template(egl::Int;name)
    @variables (fneq(t))[1:egl] (Source(t))[1:egl]

    eqs = [D(fneq) .~ Source]

    eqs = Symbolics.scalarize(reduce(vcat, Symbolics.scalarize.(eqs)))

    ODESystem(eqs,t;name)
end

function athem_excitation(fneq,Tel,laser,DOS::Spline1D,egl)
    @parameters (egrid)[1:egl] hv μ kB

    feq = FermiDirac(egrid,μ,Tel,kB)
    ftot = fneq.+feq
    ftotspl=get_interpolate(egrid,ftot)

    Δfe = fgr_electron_generation(egrid,DOS,ftotspl,hv)
    Δfh = fgr_hole_generation(egrid,DOS,ftotspl,hv)

    pc_sf = get_noparticles(Δfe,DOS,egrid) / get_noparticles(Δfh,DOS,egrid)
    Δfshape = (Δfe.*pc_sf).-Δfh
    inten = get_internalenergy_grid(Δfshape,DOS,egrid)

    return laser*Δfshape./inten
end

function fgr_hole_generation(egrid,DOS::Spline1D,ftotspl,hv)
    return DOS(egrid.+hv).*ftotspl(egrid).*(1 .-ftotspl(egrid.+hv))
end
@register_array_symbolic fgr_hole_generation(egrid::AbstractVector,DOS::Spline1D,ftotspl::Spline1D,hv::Num) begin
    size = (length(egrid),)
    eltype = eltype(egrid)
end
function fgr_electron_generation(egrid,DOS::Spline1D,ftotspl,hv)
    return DOS(egrid.-hv).*ftotspl(egrid.-hv).*(1 .-ftotspl(egrid))
end
@register_array_symbolic fgr_electron_generation(egrid::AbstractVector,DOS::Spline1D,ftotspl::Spline1D,hv::Num) begin
    size = (length(egrid),)
    eltype = eltype(egrid)
end

FermiDirac(egrid,μ::Real,Tel::Real,kB::Real) = 1 ./(exp.((egrid.-μ)./(kB*Tel)).+1)

function get_internalenergy_grid(Dis,DOS,egrid)
    integrand = Dis.*DOS(egrid).*egrid
    prob = SampledIntegralProblem(integrand,egrid)
    return solve(prob,SimpsonsRule()).u
end
@register_symbolic get_internalenergy_grid(Dis::AbstractVector,DOS::Spline1D,egrid::AbstractVector)

function get_noparticles(Dis,DOS::Spline1D,egrid)
    integrand = Dis.*DOS(egrid)
    prob = SampledIntegralProblem(integrand,egrid)
    return solve(prob,SimpsonsRule()).u
end
@register_symbolic get_noparticles(Dis::AbstractVector,DOS::Spline1D,egrid::AbstractVector)

get_interpolate(xvals::AbstractVector,yvals::AbstractVector) = Spline1D(xvals,yvals,bc="nearest")
@register_symbolic get_interpolate(xvals::AbstractVector,yvals::AbstractVector)::Spline1D

function laser()
    @parameters ϵ FWHM ϕ R
    return (0.9394372786996513(1 - R)*exp((-2.772588722239781(t^2)) / (FWHM^2))*ϕ) / (FWHM*ϵ)
end
function generate_DOS()
    return get_interpolate(collect(1:100),randn(100).*59)
end

function main()
    egl = 605
    DOS = generate_DOS()
    las = laser()
    @named test = athem_factory(DOS,las,egl)
    structural_simplify(test)
end

main()
8

I think I have somewhat traced the problem to a duplication of unknowns. When I get a method that manages to compose. The structural_simplify then fails and this seems to be because the unknown dfneq.neq appears both as a list of individual elements aka. dfneq.neq(t)[1], dfneq.neq(t)[2] etc. and as an individual dfneq.neq(t) also in the list. I am not sure what is generating this second version of the vector in the variables map.