From 117a0ef5bd2fbf441183ae51d1e5d20d932cfb2a Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 15 Mar 2022 20:58:05 -0400 Subject: [PATCH 01/19] Add allow_parameter option and make find_eq_solvables less aggressive --- src/structural_transformation/utils.jl | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl index 3d7f98a7df..348a608e6f 100644 --- a/src/structural_transformation/utils.jl +++ b/src/structural_transformation/utils.jl @@ -158,7 +158,7 @@ end ### Structural and symbolic utilities ### -function find_eq_solvables!(state::TearingState, ieq; may_be_zero=false, allow_symbolic=true) +function find_eq_solvables!(state::TearingState, ieq; may_be_zero=false, allow_symbolic=false, allow_parameter=true) fullvars = state.fullvars @unpack graph, solvable_graph = state.structure eq = equations(state)[ieq] @@ -171,7 +171,13 @@ function find_eq_solvables!(state::TearingState, ieq; may_be_zero=false, allow_s a = unwrap(a) islinear || continue if a isa Symbolic - allow_symbolic || continue + if !allow_symbolic + if allow_parameter + ModelingToolkit.isparameter(a) || continue + else + continue + end + end add_edge!(solvable_graph, ieq, j) continue end @@ -191,7 +197,7 @@ function find_eq_solvables!(state::TearingState, ieq; may_be_zero=false, allow_s end end -function find_solvables!(state::TearingState; allow_symbolic=true) +function find_solvables!(state::TearingState; allow_symbolic=false) @assert state.structure.solvable_graph === nothing eqs = equations(state) graph = state.structure.graph From 32206c984685867ce7f7c33f8ed5c03fd50a8b3d Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 15 Mar 2022 22:26:29 -0400 Subject: [PATCH 02/19] Evaluate/expand instream op at runtime --- src/systems/connectors.jl | 103 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 99 insertions(+), 4 deletions(-) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index b168326058..0f3bee18fc 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -173,14 +173,39 @@ function collect_instream!(set, expr, occurs=false) return occurs end -positivemax(m, ::Any; tol=nothing)= max(m, something(tol, 1e-8)) +#positivemax(m, ::Any; tol=nothing)= max(m, something(tol, 1e-8)) +#_positivemax(m, tol) = ifelse((-tol <= m) & (m <= tol), ((3 * tol - m) * (tol + m)^3)/(16 * tol^3) + tol, max(m, tol)) +function _positivemax(m, si) + T = typeof(m) + relativeTolerance = 1e-4 + nominal = one(T) + eps = relativeTolerance * nominal + alpha = if si > eps + one(T) + else + if si > 0 + (si/eps)^2*(3-2* si/eps) + else + zero(T) + end + end + alpha * max(m, 0) + (1-alpha)*eps +end +@register _positivemax(m, tol) +positivemax(m, ::Any; tol=nothing) = _positivemax(m, tol) +mydiv(num, den) = if den == 0 + error() +else + num / den +end +@register mydiv(n, d) function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10) subsys = get_systems(sys) isempty(subsys) && return sys # post order traversal - @set! sys.systems = map(s->expand_connections(s, debug=debug), subsys) + @set! sys.systems = map(s->expand_connections(s, debug=debug, tol=tol), subsys) outer_connectors = Symbol[] for s in subsys @@ -288,6 +313,67 @@ function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10) return sys end +@generated function _instream_split(::Val{inner_n}, ::Val{outer_n}, vars::NTuple{N}) where {inner_n, outer_n, N} + #instream_rt(innerfvs..., innersvs..., outerfvs..., outersvs...) + ret = Expr(:tuple) + # mj.c.m_flow + inner_f = :(Base.@ntuple $inner_n i -> vars[i]) + offset = inner_n + inner_s = :(Base.@ntuple $inner_n i -> vars[$offset+i]) + offset += inner_n + # ck.m_flow + outer_f = :(Base.@ntuple $outer_n i -> vars[$offset+i]) + offset += outer_n + outer_s = :(Base.@ntuple $outer_n i -> vars[$offset+i]) + Expr(:tuple, inner_f, inner_s, outer_f, outer_s) +end + +function instream_rt(ins::Val{inner_n}, outs::Val{outer_n}, vars::Vararg{Any,N}) where {inner_n, outer_n, N} + @assert N == 2*(inner_n + outer_n) + + # inner: mj.c.m_flow + # outer: ck.m_flow + inner_f, inner_s, outer_f, outer_s = _instream_split(ins, outs, vars) + + T = float(first(inner_f)) + si = zero(T) + num = den = zero(T) + for f in inner_f + si += max(-f, 0) + end + for f in outer_f + si += max(f, 0) + end + #for (f, s) in zip(inner_f, inner_s) + for j in 1:inner_n + @inbounds f = inner_f[j] + @inbounds s = inner_s[j] + num += _positivemax(-f, si) * s + den += _positivemax(-f, si) + end + #for (f, s) in zip(outer_f, outer_s) + for j in 1:outer_n + @inbounds f = outer_f[j] + @inbounds s = outer_s[j] + num += _positivemax(-f, si) * s + den += _positivemax(-f, si) + end + return num / den + #= + si = sum(max(-mj.c.m_flow,0) for j in cat(1,1:i-1, i+1:N)) + + sum(max(ck.m_flow ,0) for k in 1:M) + + inStream(mi.c.h_outflow) = + (sum(positiveMax(-mj.c.m_flow,si)*mj.c.h_outflow) + + sum(positiveMax(ck.m_flow,s_i)*inStream(ck.h_outflow)))/ + (sum(positiveMax(-mj.c.m_flow,s_i)) + + sum(positiveMax(ck.m_flow,s_i))) + for j in 1:N and i <> j and mj.c.m_flow.min < 0, + for k in 1:M and ck.m_flow.max > 0 + =# +end +SymbolicUtils.promote_symtype(::typeof(instream_rt), ::Vararg) = Real + function expand_instream(instream_eqs, instream_exprs, connects; debug=false, tol) sub = Dict() seen = Set() @@ -339,6 +425,15 @@ function expand_instream(instream_eqs, instream_exprs, connects; debug=false, to fv = flowvar(first(connectors)) i = findfirst(c->nameof(c) === connector_name, inner_sc) if i !== nothing + # mj.c.m_flow + innerfvs = [unwrap(states(s, fv)) for (j, s) in enumerate(inner_sc) if j != i] + innersvs = [unwrap(states(s, sv)) for (j, s) in enumerate(inner_sc) if j != i] + # ck.m_flow + outerfvs = [unwrap(states(s, fv)) for s in outer_sc] + outersvs = [instream(states(s, fv)) for s in outer_sc] + + sub[ex] = term(instream_rt, Val(length(innerfvs)), Val(length(outerfvs)), innerfvs..., innersvs..., outerfvs..., outersvs...) + #= si = isempty(outer_sc) ? 0 : sum(s->max(states(s, fv), 0), outer_sc) for j in 1:n_inners; j == i && continue f = states(inner_sc[j], fv) @@ -359,9 +454,9 @@ function expand_instream(instream_eqs, instream_exprs, connects; debug=false, to den += tmp num += tmp * instream(states(outer_sc[k], sv)) end - sub[ex] = num / den + sub[ex] = mydiv(num, den) + =# end - end end From e47fb185f4bd1a32782bec22aebe965887f6b7b4 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Fri, 25 Mar 2022 15:28:37 -0400 Subject: [PATCH 03/19] Ad Hoc "fix" --- src/systems/connectors.jl | 41 ++++++++++++++++++++++++++++++++------- 1 file changed, 34 insertions(+), 7 deletions(-) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index 0f3bee18fc..69d368fdab 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -56,7 +56,18 @@ end Connection(syss) = Connection(inners=syss) get_systems(c::Connection) = c.inners function Base.in(e::Symbol, c::Connection) - any(k->nameof(k) === e, c.inners) || any(k->nameof(k) === e, c.outers) + (c.inners !== nothing && any(k->nameof(k) === e, c.inners)) || + (c.outers !== nothing && any(k->nameof(k) === e, c.outers)) +end + +function renamespace(sym::Symbol, connection::Connection) + inners = connection.inners === nothing ? [] : renamespace.(sym, connection.inners) + if connection.outers !== nothing + for o in connection.outers + push!(inners, renamespace(sym, o)) + end + end + Connection(;inners=inners) end const EMPTY_VEC = [] @@ -200,12 +211,14 @@ else end @register mydiv(n, d) -function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10) +function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10, + rename=Ref{Union{Nothing,Tuple{Symbol,Int}}}(nothing), stream_connects=[]) subsys = get_systems(sys) isempty(subsys) && return sys # post order traversal - @set! sys.systems = map(s->expand_connections(s, debug=debug, tol=tol), subsys) + @set! sys.systems = map(s->expand_connections(s, debug=debug, tol=tol, + rename=rename, stream_connects=stream_connects), subsys) outer_connectors = Symbol[] for s in subsys @@ -306,7 +319,21 @@ function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10) end # stream variables - stream_connects = filter(isstreamconnection, narg_connects) + if rename[] !== nothing + name, depth = rename[] + nc = length(stream_connects) + for i in nc-depth+1:nc + stream_connects[i] = renamespace(:room, stream_connects[i]) + end + end + nsc = 0 + for c in narg_connects + if isstreamconnection(c) + push!(stream_connects, c) + nsc += 1 + end + end + rename[] = nameof(sys), nsc instream_eqs, additional_eqs = expand_instream(instream_eqs, instream_exprs, stream_connects; debug=debug, tol=tol) @set! sys.eqs = [eqs; instream_eqs; additional_eqs] @@ -388,8 +415,8 @@ function expand_instream(instream_eqs, instream_exprs, connects; debug=false, to connectors = Iterators.flatten((connect.inners, connect.outers)) # stream variable sv = getproperty(first(connectors), streamvar_name; namespace=false) - inner_sc = connect.inners - outer_sc = connect.outers + inner_sc = something(connect.inners, EMPTY_VEC) + outer_sc = something(connect.outers, EMPTY_VEC) n_outers = length(outer_sc) n_inners = length(inner_sc) @@ -463,7 +490,7 @@ function expand_instream(instream_eqs, instream_exprs, connects; debug=false, to # additional equations additional_eqs = Equation[] for c in connects - outer_sc = c.outers + outer_sc = something(c.outers, EMPTY_VEC) isempty(outer_sc) && continue inner_sc = c.inners n_outers = length(outer_sc) From 417530b57d18bbc9b105ae21e6bbe393d991573a Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Mon, 4 Apr 2022 16:01:42 -0400 Subject: [PATCH 04/19] WIP generate_connection_set --- src/systems/connectors.jl | 60 ++++++++++++++++++++++++++++++++------- 1 file changed, 49 insertions(+), 11 deletions(-) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index 69d368fdab..b1c1e869b8 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -211,21 +211,13 @@ else end @register mydiv(n, d) -function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10, - rename=Ref{Union{Nothing,Tuple{Symbol,Int}}}(nothing), stream_connects=[]) - subsys = get_systems(sys) - isempty(subsys) && return sys - - # post order traversal - @set! sys.systems = map(s->expand_connections(s, debug=debug, tol=tol, - rename=rename, stream_connects=stream_connects), subsys) - +function generate_isouter(sys::AbstractSystem) outer_connectors = Symbol[] - for s in subsys + for s in get_systems(sys) n = nameof(s) isconnector(s) && push!(outer_connectors, n) end - isouter = let outer_connectors=outer_connectors + let outer_connectors=outer_connectors function isouter(sys)::Bool s = string(nameof(sys)) isconnector(sys) || error("$s is not a connector!") @@ -234,7 +226,53 @@ function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10, parent_name in outer_connectors end end +end + +struct ConnectionSet + set::Vector{Pair{Any, Bool}} # var => isouter +end + +function generate_connection_set(sys::AbstractSystem, namespace=nothing) + subsys = get_systems(sys) + isempty(subsys) && return sys + + isouter = generate_isouter(sys) + eqs′ = get_eqs(sys) + eqs = Equation[] + #instream_eqs = Equation[] + #instream_exprs = [] + cts = [] # connections + for eq in eqs′ + if eq.lhs isa Connection + push!(cts, get_systems(eq.rhs)) + #elseif collect_instream!(instream_exprs, eq) + # push!(instream_eqs, eq) + else + push!(eqs, eq) # split connections and equations + end + end + + # if there are no connections, we are done + isempty(cts) && return sys + + set = Pair{Any, Bool}[] + + + # pre order traversal + namespace = renamespace(nameof(sys), namespace) + @set! sys.systems = map(Base.Fix2(generate_connection_set, namespace), subsys) +end + +function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10, + rename=Ref{Union{Nothing,Tuple{Symbol,Int}}}(nothing), stream_connects=[]) + subsys = get_systems(sys) + isempty(subsys) && return sys + + # post order traversal + @set! sys.systems = map(s->expand_connections(s, debug=debug, tol=tol, + rename=rename, stream_connects=stream_connects), subsys) + isouter = generate_isouter(sys) sys = flatten(sys) eqs′ = get_eqs(sys) eqs = Equation[] From 74abd4707981da9d6093fc2f1ae558c1e6258ff9 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Mon, 4 Apr 2022 16:02:00 -0400 Subject: [PATCH 05/19] Add Modelica's example in test --- test/components.jl | 47 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/test/components.jl b/test/components.jl index 188ef55cef..e0311aca5f 100644 --- a/test/components.jl +++ b/test/components.jl @@ -138,3 +138,50 @@ end @unpack foo = goo @test ModelingToolkit.defaults(goo)[foo.a] == 3 @test ModelingToolkit.defaults(goo)[foo.b] == 300 + +#= +model Circuit + Ground ground; + Load load; + Resistor resistor; +equation + connect(load.p , ground.p); + connect(resistor.p, ground.p); +end Circuit; +model Load + extends TwoPin; + Resistor resistor; +equation + connect(p, resistor.p); + connect(resistor.n, n); +end Load; +=# + +function Load(;name) + R = 1 + @named p = Pin() + @named n = Pin() + @named resistor = Resistor(R=R) + eqs = [ + connect(p, resistor.p); + connect(resistor.n, n); + ] + @named sys = ODESystem(eqs, t) + compose(sys, [p, n, resistor]; name=name) +end + +function Circuit(;name) + R = 1 + @named ground = Ground() + @named load = Load() + @named resistor = Resistor(R=R) + eqs = [ + connect(load.p , ground.g); + connect(resistor.p, ground.g); + ] + @named sys = ODESystem(eqs, t) + compose(sys, [ground, resistor, load]; name=name) +end + +@named foo = Circuit() +@test_broken structural_simplify(foo) isa AbstractSystem From acadb7eca4bd3f45bb0471f1788880803a818737 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Mon, 4 Apr 2022 17:53:20 -0400 Subject: [PATCH 06/19] Add generate_connection_set that behaves like Modelica --- src/systems/abstractsystem.jl | 1 + src/systems/connectors.jl | 106 ++++++++++++++++++++++++++++++++-- 2 files changed, 102 insertions(+), 5 deletions(-) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 079111ec89..519c3c57d4 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -380,6 +380,7 @@ renamespace(sys, eq::Equation) = namespace_equation(eq, sys) renamespace(names::AbstractVector, x) = foldr(renamespace, names, init=x) function renamespace(sys, x) + sys === nothing && return x x = unwrap(x) if x isa Symbolic let scope = getmetadata(x, SymScope, LocalScope()) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index b1c1e869b8..a52a5b1a1c 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -232,8 +232,43 @@ struct ConnectionSet set::Vector{Pair{Any, Bool}} # var => isouter end -function generate_connection_set(sys::AbstractSystem, namespace=nothing) +function Base.show(io::IO, c::ConnectionSet) + print(io, "<") + for i in 1:length(c.set)-1 + v, isouter = c.set[i] + print(io, v, "::", isouter ? "outer" : "inner", ", ") + end + v, isouter = last(c.set) + print(io, v, "::", isouter ? "outer" : "inner", ">") +end + +@noinline connection_error(ss) = error("Different types of connectors are in one conenction statement: <$(map(nameof, ss))>") + +function connection2set!(connectionsets, namespace, ss, isouter) + sts1 = Set(states(first(ss))) + T = Pair{Any,Bool} + csets = [T[] for _ in 1:length(ss)] + for (i, s) in enumerate(ss) + sts = states(s) + i != 1 && ((length(sts1) == length(sts) && all(Base.Fix2(in, sts1), sts)) || connection_error(ss)) + io = isouter(s) + for (j, v) in enumerate(sts) + push!(csets[j], T(states(renamespace(namespace, s), v), io)) + end + end + for cset in csets + vtype = get_connection_type(first(cset)[1]) + for k in 2:length(cset) + vtype === get_connection_type(cset[k][1]) || connection_error(ss) + end + push!(connectionsets, ConnectionSet(cset)) + end +end + +generate_connection_set(sys::AbstractSystem) = (connectionsets = ConnectionSet[]; generate_connection_set!(connectionsets, sys::AbstractSystem); connectionsets) +function generate_connection_set!(connectionsets, sys::AbstractSystem, namespace=nothing) subsys = get_systems(sys) + # no connectors if there are no subsystems isempty(subsys) && return sys isouter = generate_isouter(sys) @@ -252,16 +287,77 @@ function generate_connection_set(sys::AbstractSystem, namespace=nothing) end end + if namespace !== nothing + # Except for the top level, all connectors are eventually inside + # connectors. + T = Pair{Any,Bool} + for s in subsys + isconnector(s) || continue + for v in states(s) + Flow === get_connection_type(v) || continue + push!(connectionsets, ConnectionSet([T(renamespace(namespace, states(s, v)), false)])) + end + end + end + # if there are no connections, we are done isempty(cts) && return sys - set = Pair{Any, Bool}[] - + for ct in cts + connection2set!(connectionsets, namespace, ct, isouter) + end # pre order traversal - namespace = renamespace(nameof(sys), namespace) - @set! sys.systems = map(Base.Fix2(generate_connection_set, namespace), subsys) + @set! sys.systems = map(s->generate_connection_set!(connectionsets, s, renamespace(namespace, nameof(s))), subsys) end +#= + + +{} + + +{} + + +{} + + +{} + + +{} + + +{} + + +{} + + + +{, } + + +{, } + + +{, } + + +{, } + + +{, } + + +{, } + + +{, } + + +{, } +=# function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10, rename=Ref{Union{Nothing,Tuple{Symbol,Int}}}(nothing), stream_connects=[]) From ebbcfeef094163b10c643aa45b02fdbedfcbd958 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Mon, 4 Apr 2022 18:19:12 -0400 Subject: [PATCH 07/19] Add Modelica compatible generate_connection_equations --- src/systems/connectors.jl | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index a52a5b1a1c..2e8954f072 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -359,6 +359,42 @@ end {, } =# +function Base.merge(csets::AbstractVector{<:ConnectionSet}) + mcsets = ConnectionSet[] + # FIXME: this is O(m n^3) + for cset in csets + idx = findfirst(mcset->any(s->any(z->isequal(z, s), cset.set), mcset.set), mcsets) + if idx === nothing + push!(mcsets, cset) + else + union!(mcsets[idx].set, cset.set) + end + end + mcsets +end + +function generate_connection_equations(csets::AbstractVector{<:ConnectionSet}) + eqs = Equation[] + for cset in csets + vtype = get_connection_type(cset.set[1][1]) + vtype === Stream && continue + if vtype === Flow + rhs = 0 + for (v, isouter) in cset.set + rhs += isouter ? -v : v + end + push!(eqs, 0 ~ rhs) + else # Equality + base = cset.set[1][1] + for i in 2:length(cset.set) + v = cset.set[i][1] + push!(eqs, base ~ v) + end + end + end + eqs +end + function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10, rename=Ref{Union{Nothing,Tuple{Symbol,Int}}}(nothing), stream_connects=[]) subsys = get_systems(sys) From 0b0f7844f9d83d6752b7f99ea3116cc4143ff153 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 5 Apr 2022 19:17:59 -0400 Subject: [PATCH 08/19] Clean up --- src/systems/connectors.jl | 121 ++++++++++++++++---------------------- 1 file changed, 50 insertions(+), 71 deletions(-) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index 2e8954f072..f57e61e253 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -228,44 +228,66 @@ function generate_isouter(sys::AbstractSystem) end end +struct LazyNamespace + namespace::Union{Nothing,Symbol} + sys +end + +Base.copy(l::LazyNamespace) = renamespace(l.namespace, l.sys) +Base.nameof(l::LazyNamespace) = renamespace(l.namespace, nameof(l.sys)) + +struct ConnectionElement + sys::LazyNamespace + v + isouter::Bool +end +Base.hash(l::ConnectionElement, salt::UInt) = hash(nameof(l.sys)) ⊻ hash(l.v) ⊻ hash(l.isouter) ⊻ salt +Base.isequal(l1::ConnectionElement, l2::ConnectionElement) = l1 == l2 +Base.:(==)(l1::ConnectionElement, l2::ConnectionElement) = nameof(l1.sys) == nameof(l2.sys) && isequal(l1.v, l2.v) && l1.isouter == l2.isouter +function namespaced_var(l::ConnectionElement) + @unpack sys, v = l + states(copy(sys), v) +end + struct ConnectionSet - set::Vector{Pair{Any, Bool}} # var => isouter + set::Vector{ConnectionElement} # namespace.sys, var, isouter end function Base.show(io::IO, c::ConnectionSet) print(io, "<") for i in 1:length(c.set)-1 - v, isouter = c.set[i] - print(io, v, "::", isouter ? "outer" : "inner", ", ") + @unpack sys, v, isouter = c.set[i] + print(io, nameof(sys), ".", v, "::", isouter ? "outer" : "inner", ", ") end - v, isouter = last(c.set) - print(io, v, "::", isouter ? "outer" : "inner", ">") + @unpack sys, v, isouter = last(c.set) + print(io, nameof(sys), ".", v, "::", isouter ? "outer" : "inner", ">") end @noinline connection_error(ss) = error("Different types of connectors are in one conenction statement: <$(map(nameof, ss))>") function connection2set!(connectionsets, namespace, ss, isouter) + nn = map(nameof, ss) sts1 = Set(states(first(ss))) - T = Pair{Any,Bool} - csets = [T[] for _ in 1:length(ss)] + T = ConnectionElement + csets = [T[] for _ in 1:length(sts1)] for (i, s) in enumerate(ss) sts = states(s) i != 1 && ((length(sts1) == length(sts) && all(Base.Fix2(in, sts1), sts)) || connection_error(ss)) io = isouter(s) for (j, v) in enumerate(sts) - push!(csets[j], T(states(renamespace(namespace, s), v), io)) + push!(csets[j], T(LazyNamespace(namespace, s), v, io)) end end for cset in csets - vtype = get_connection_type(first(cset)[1]) + vtype = get_connection_type(first(cset).v) for k in 2:length(cset) - vtype === get_connection_type(cset[k][1]) || connection_error(ss) + vtype === get_connection_type(cset[k].v) || connection_error(ss) end push!(connectionsets, ConnectionSet(cset)) end end -generate_connection_set(sys::AbstractSystem) = (connectionsets = ConnectionSet[]; generate_connection_set!(connectionsets, sys::AbstractSystem); connectionsets) +generate_connection_set(sys::AbstractSystem) = (connectionsets = ConnectionSet[]; (generate_connection_set!(connectionsets, sys::AbstractSystem), connectionsets)) function generate_connection_set!(connectionsets, sys::AbstractSystem, namespace=nothing) subsys = get_systems(sys) # no connectors if there are no subsystems @@ -290,12 +312,12 @@ function generate_connection_set!(connectionsets, sys::AbstractSystem, namespace if namespace !== nothing # Except for the top level, all connectors are eventually inside # connectors. - T = Pair{Any,Bool} + T = ConnectionElement for s in subsys isconnector(s) || continue for v in states(s) Flow === get_connection_type(v) || continue - push!(connectionsets, ConnectionSet([T(renamespace(namespace, states(s, v)), false)])) + push!(connectionsets, ConnectionSet([T(LazyNamespace(namespace, s), v, false)])) end end end @@ -309,61 +331,14 @@ function generate_connection_set!(connectionsets, sys::AbstractSystem, namespace # pre order traversal @set! sys.systems = map(s->generate_connection_set!(connectionsets, s, renamespace(namespace, nameof(s))), subsys) + @set! sys.eqs = eqs end -#= - - -{} - - -{} - - -{} - - -{} - - -{} - - -{} - - -{} - - - -{, } - - -{, } - - -{, } - - -{, } - - -{, } - - -{, } - - -{, } - - -{, } -=# function Base.merge(csets::AbstractVector{<:ConnectionSet}) mcsets = ConnectionSet[] # FIXME: this is O(m n^3) for cset in csets - idx = findfirst(mcset->any(s->any(z->isequal(z, s), cset.set), mcset.set), mcsets) + idx = findfirst(mcset->any(s->any(z->z == s, cset.set), mcset.set), mcsets) if idx === nothing push!(mcsets, cset) else @@ -373,26 +348,30 @@ function Base.merge(csets::AbstractVector{<:ConnectionSet}) mcsets end -function generate_connection_equations(csets::AbstractVector{<:ConnectionSet}) +function generate_connection_equations_and_stream_connections(csets::AbstractVector{<:ConnectionSet}) eqs = Equation[] + stream_connections = ConnectionSet[] for cset in csets - vtype = get_connection_type(cset.set[1][1]) - vtype === Stream && continue - if vtype === Flow + vtype = get_connection_type(cset.set[1].v) + if vtype === Stream + push!(stream_connections, cset) + continue + elseif vtype === Flow rhs = 0 - for (v, isouter) in cset.set - rhs += isouter ? -v : v + for ele in cset.set + v = namespaced_var(ele) + rhs += ele.isouter ? -v : v end push!(eqs, 0 ~ rhs) else # Equality - base = cset.set[1][1] + base = namespaced_var(cset.set[1]) for i in 2:length(cset.set) - v = cset.set[i][1] + v = namespaced_var(cset.set[i]) push!(eqs, base ~ v) end end end - eqs + eqs, stream_connections end function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10, From 73c39b27a7cb98a8d9f024e968e46ee3e1c5b5eb Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 5 Apr 2022 21:27:20 -0400 Subject: [PATCH 09/19] Add expand_instream2 --- src/systems/connectors.jl | 160 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 155 insertions(+), 5 deletions(-) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index f57e61e253..ef17d46c10 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -244,10 +244,8 @@ end Base.hash(l::ConnectionElement, salt::UInt) = hash(nameof(l.sys)) ⊻ hash(l.v) ⊻ hash(l.isouter) ⊻ salt Base.isequal(l1::ConnectionElement, l2::ConnectionElement) = l1 == l2 Base.:(==)(l1::ConnectionElement, l2::ConnectionElement) = nameof(l1.sys) == nameof(l2.sys) && isequal(l1.v, l2.v) && l1.isouter == l2.isouter -function namespaced_var(l::ConnectionElement) - @unpack sys, v = l - states(copy(sys), v) -end +namespaced_var(l::ConnectionElement) = states(l, l.v) +states(l::ConnectionElement, v) = states(copy(l.sys), v) struct ConnectionSet set::Vector{ConnectionElement} # namespace.sys, var, isouter @@ -550,9 +548,161 @@ function instream_rt(ins::Val{inner_n}, outs::Val{outer_n}, vars::Vararg{Any,N}) end SymbolicUtils.promote_symtype(::typeof(instream_rt), ::Vararg) = Real +function expand_instream2(csets::AbstractVector{<:ConnectionSet}, sys::AbstractSystem, namespace=nothing, debug=false, tol=1e-8) + subsys = get_systems(sys) + # no connectors if there are no subsystems + isempty(subsys) && return sys + # post order traversal + @set! sys.systems = map(s->expand_instream2(csets, s, renamespace(namespace, nameof(s)), debug, tol), subsys) + + eqs′ = get_eqs(sys) + eqs = Equation[] + instream_eqs = Equation[] + instream_exprs = [] + for eq in eqs′ + if collect_instream!(instream_exprs, eq) + push!(instream_eqs, eq) + else + push!(eqs, eq) # split connections and equations + end + end + + + sub = Dict() + for ex in instream_exprs + sv = only(arguments(ex)) + full_name_sv = renamespace(namespace, sv) + + #cidx = findfirst(c->any(v->, ), ) + cidx = -1 + idx_in_set = -1 + for (i, c) in enumerate(csets) + for (j, v) in enumerate(c.set) + if isequal(namespaced_var(v), full_name_sv) + cidx = i + idx_in_set = j + end + end + end + cidx < 0 && error("$sv is not a variable inside stream connectors") + cset = csets[cidx].set + + connectors = Vector{Any}(undef, length(cset)) + n_inners = n_outers = 0 + for (i, e) in enumerate(cset) + connectors[i] = e.sys.sys + if e.isouter + n_outers += 1 + else + n_inners += 1 + end + end + if n_inners == 1 && n_outers == 0 + sub[ex] = sv + elseif n_inners == 2 && n_outers == 0 + other = idx_in_set == 1 ? 2 : 1 + sub[ex] = states(cset[other], sv) + elseif n_inners == 1 && n_outers == 1 + if !cset[idx_in_set].isouter + other = idx_in_set == 1 ? 2 : 1 + outerstream = states(cset[other], sv) + sub[ex] = outerstream + end + else + if !cset[idx_in_set].isouter + fv = flowvar(first(connectors)) + # mj.c.m_flow + innerfvs = [unwrap(states(s, fv)) for (j, s) in enumerate(cset) if j != idx_in_set && !s.isouter] + innersvs = [unwrap(states(s, sv)) for (j, s) in enumerate(cset) if j != idx_in_set && !s.isouter] + # ck.m_flow + outerfvs = [unwrap(states(s, fv)) for s in cset if s.isouter] + outersvs = [unwrap(states(s, sv)) for s in cset if s.isouter] + + sub[ex] = term(instream_rt, Val(length(innerfvs)), Val(length(outerfvs)), innerfvs..., innersvs..., outerfvs..., outersvs...) + end + end + end + + # additional equations + additional_eqs = Equation[] + csets = filter(cset->any(e->e.sys.namespace === namespace, cset.set), csets) + for cset′ in csets + cset = cset′.set + connectors = Vector{Any}(undef, length(cset)) + n_inners = n_outers = 0 + for (i, e) in enumerate(cset) + connectors[i] = e.sys.sys + if e.isouter + n_outers += 1 + else + n_inners += 1 + end + end + iszero(n_outers) && continue + connector_representative = first(cset).sys.sys + fv = flowvar(connector_representative) + sv = first(cset).v + vtype = get_connection_type(sv) + vtype === Stream || continue + if n_inners == 1 && n_outers == 1 + push!(additional_eqs, states(cset[1], sv) ~ states(cset[2], sv)) + elseif n_inners == 0 && n_outers == 2 + # we don't expand `instream` in this case. + v1 = states(cset[1], sv) + v2 = states(cset[2], sv) + push!(additional_eqs, v1 ~ instream(v2)) + push!(additional_eqs, v2 ~ instream(v1)) + else + sq = 0 + s_inners = (s for s in cset if !s.isouter) + s_outers = (s for s in cset if s.isouter) + for (q, oscq) in enumerate(s_outers) + sq += sum(s->max(-states(s, fv), 0), s_inners) + for (k, s) in enumerate(s_outers); k == q && continue + f = states(s, fv) + sq += max(f, 0) + end + + num = 0 + den = 0 + for s in s_inners + f = states(s, fv) + tmp = positivemax(-f, sq; tol=tol) + den += tmp + num += tmp * states(s, sv) + end + for (k, s) in enumerate(s_outers); k == q && continue + f = states(s, fv) + tmp = positivemax(f, sq; tol=tol) + den += tmp + num += tmp * instream(states(s, sv)) + end + push!(additional_eqs, states(oscq, sv) ~ num / den) + end + end + end + + instream_eqs = map(Base.Fix2(substitute, sub), instream_eqs) + if debug + println("===========BEGIN=============") + println("Expanded equations:") + for eq in instream_eqs + print_with_indent(4, eq) + end + if !isempty(additional_eqs) + println("Additional equations:") + for eq in additional_eqs + print_with_indent(4, eq) + end + end + println("============END==============") + end + + @set! sys.eqs = [eqs; instream_eqs; additional_eqs] +end + function expand_instream(instream_eqs, instream_exprs, connects; debug=false, tol) sub = Dict() - seen = Set() for ex in instream_exprs var = only(arguments(ex)) connector_name, streamvar_name = split_sys_var(var) From 74426792522bf55414368d2427fae34cf65a0775 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Thu, 7 Apr 2022 18:15:41 -0400 Subject: [PATCH 10/19] Fix instream typo --- src/systems/connectors.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index ef17d46c10..67194c1cd0 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -487,7 +487,7 @@ function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10, return sys end -@generated function _instream_split(::Val{inner_n}, ::Val{outer_n}, vars::NTuple{N}) where {inner_n, outer_n, N} +@generated function _instream_split(::Val{inner_n}, ::Val{outer_n}, vars::NTuple{N,Any}) where {inner_n, outer_n, N} #instream_rt(innerfvs..., innersvs..., outerfvs..., outersvs...) ret = Expr(:tuple) # mj.c.m_flow @@ -756,7 +756,7 @@ function expand_instream(instream_eqs, instream_exprs, connects; debug=false, to innersvs = [unwrap(states(s, sv)) for (j, s) in enumerate(inner_sc) if j != i] # ck.m_flow outerfvs = [unwrap(states(s, fv)) for s in outer_sc] - outersvs = [instream(states(s, fv)) for s in outer_sc] + outersvs = [instream(states(s, sv)) for s in outer_sc] sub[ex] = term(instream_rt, Val(length(innerfvs)), Val(length(outerfvs)), innerfvs..., innersvs..., outerfvs..., outersvs...) #= From e767cfb9d9bd76f89089213829d5048ac81dcf8a Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Thu, 7 Apr 2022 19:06:30 -0400 Subject: [PATCH 11/19] Use the new connection set generation in expand_connections --- src/systems/connectors.jl | 39 ++++++++++++++------------------------- 1 file changed, 14 insertions(+), 25 deletions(-) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index 67194c1cd0..d593d5ac41 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -285,7 +285,7 @@ function connection2set!(connectionsets, namespace, ss, isouter) end end -generate_connection_set(sys::AbstractSystem) = (connectionsets = ConnectionSet[]; (generate_connection_set!(connectionsets, sys::AbstractSystem), connectionsets)) +generate_connection_set(sys::AbstractSystem) = (connectionsets = ConnectionSet[]; (generate_connection_set!(connectionsets, sys::AbstractSystem), merge(connectionsets))) function generate_connection_set!(connectionsets, sys::AbstractSystem, namespace=nothing) subsys = get_systems(sys) # no connectors if there are no subsystems @@ -293,7 +293,8 @@ function generate_connection_set!(connectionsets, sys::AbstractSystem, namespace isouter = generate_isouter(sys) eqs′ = get_eqs(sys) - eqs = Equation[] + #eqs = Equation[] + #instream_eqs = Equation[] #instream_exprs = [] cts = [] # connections @@ -303,7 +304,7 @@ function generate_connection_set!(connectionsets, sys::AbstractSystem, namespace #elseif collect_instream!(instream_exprs, eq) # push!(instream_eqs, eq) else - push!(eqs, eq) # split connections and equations + #push!(eqs, eq) # split connections and equations end end @@ -329,7 +330,7 @@ function generate_connection_set!(connectionsets, sys::AbstractSystem, namespace # pre order traversal @set! sys.systems = map(s->generate_connection_set!(connectionsets, s, renamespace(namespace, nameof(s))), subsys) - @set! sys.eqs = eqs + #@set! sys.eqs = eqs end function Base.merge(csets::AbstractVector{<:ConnectionSet}) @@ -372,13 +373,20 @@ function generate_connection_equations_and_stream_connections(csets::AbstractVec eqs, stream_connections end -function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10, +function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10) + sys, csets = generate_connection_set(sys) + ceqs, _ = generate_connection_equations_and_stream_connections(csets) + sys = _expand_connections(sys; debug=debug, tol=tol) + @set! sys.eqs = [equations(sys); ceqs] +end + +function _expand_connections(sys::AbstractSystem; debug=false, tol=1e-10, rename=Ref{Union{Nothing,Tuple{Symbol,Int}}}(nothing), stream_connects=[]) subsys = get_systems(sys) isempty(subsys) && return sys # post order traversal - @set! sys.systems = map(s->expand_connections(s, debug=debug, tol=tol, + @set! sys.systems = map(s->_expand_connections(s, debug=debug, tol=tol, rename=rename, stream_connects=stream_connects), subsys) isouter = generate_isouter(sys) @@ -446,25 +454,6 @@ function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10, length(dups) == 0 || error("Connection([$(join(map(nameof, allconnectors), ", "))]) has duplicated connections: [$(join(collect(dups), ", "))].") end - if debug - println("============BEGIN================") - println("Connections for [$(nameof(sys))]:") - foreach(Base.Fix1(print_with_indent, 4), narg_connects) - end - - connection_eqs = Equation[] - for c in narg_connects - ceqs = connect(c) - debug && append!(connection_eqs, ceqs) - append!(eqs, ceqs) - end - - if debug - println("Connection equations:") - foreach(Base.Fix1(print_with_indent, 4), connection_eqs) - println("=============END=================") - end - # stream variables if rename[] !== nothing name, depth = rename[] From 90ddc86802f4f493f5251797e6f4952fee5920ef Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Thu, 7 Apr 2022 19:18:33 -0400 Subject: [PATCH 12/19] Fix connect equations generation for array valued connectors --- src/systems/connectors.jl | 7 ++++++- test/stream_connectors.jl | 16 ++++++++-------- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index d593d5ac41..fa1b94764d 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -350,8 +350,13 @@ end function generate_connection_equations_and_stream_connections(csets::AbstractVector{<:ConnectionSet}) eqs = Equation[] stream_connections = ConnectionSet[] + for cset in csets - vtype = get_connection_type(cset.set[1].v) + v = cset.set[1].v + if hasmetadata(v, Symbolics.GetindexParent) + v = getparent(v) + end + vtype = get_connection_type(v) if vtype === Stream push!(stream_connections, cset) continue diff --git a/test/stream_connectors.jl b/test/stream_connectors.jl index 39b6368aa8..fccd934ec3 100644 --- a/test/stream_connectors.jl +++ b/test/stream_connectors.jl @@ -228,11 +228,11 @@ end @named simple = ODESystem([connect(vp1, vp2, vp3)], t) sys = expand_connections(compose(simple, [vp1, vp2, vp3])) -@test equations(sys) == [ - vp1.v[1] ~ vp2.v[1] - vp1.v[2] ~ vp2.v[2] - vp1.v[1] ~ vp3.v[1] - vp1.v[2] ~ vp3.v[2] - 0 ~ -vp1.i[1] - vp2.i[1] - vp3.i[1] - 0 ~ -vp1.i[2] - vp2.i[2] - vp3.i[2] - ] +@test sort(equations(sys), by=string) == sort([ + vp1.v[1] ~ vp2.v[1] + vp1.v[2] ~ vp2.v[2] + vp1.v[1] ~ vp3.v[1] + vp1.v[2] ~ vp3.v[2] + 0 ~ -vp1.i[1] - vp2.i[1] - vp3.i[1] + 0 ~ -vp1.i[2] - vp2.i[2] - vp3.i[2] + ], by=string) From 58f6dd9ca12706ea9548a8bd6b528cc1c77fcb28 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Fri, 8 Apr 2022 09:34:35 -0400 Subject: [PATCH 13/19] Relax constrains of compose --- src/systems/abstractsystem.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 519c3c57d4..0253866d5b 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -1035,14 +1035,14 @@ $(SIGNATURES) compose multiple systems together. The resulting system would inherit the first system's name. """ -function compose(sys::AbstractSystem, systems::AbstractArray{<:AbstractSystem}; name=nameof(sys)) +function compose(sys::AbstractSystem, systems::AbstractArray; name=nameof(sys)) nsys = length(systems) - nsys >= 1 || throw(ArgumentError("There must be at least 1 subsystem. Got $nsys subsystems.")) + nsys == 0 && return sys @set! sys.name = name @set! sys.systems = [get_systems(sys); systems] return sys end -compose(syss::AbstractSystem...; name=nameof(first(syss))) = compose(first(syss), collect(syss[2:end]); name=name) +compose(syss...; name=nameof(first(syss))) = compose(first(syss), collect(syss[2:end]); name=name) Base.:(∘)(sys1::AbstractSystem, sys2::AbstractSystem) = compose(sys1, sys2) UnPack.unpack(sys::ModelingToolkit.AbstractSystem, ::Val{p}) where p = getproperty(sys, p; namespace=false) From f179992bc403480a86f35abb280b354da4013592 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Fri, 8 Apr 2022 20:55:06 -0400 Subject: [PATCH 14/19] Fix invalidate_cache --- src/systems/abstractsystem.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 0253866d5b..0dda0df4a2 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -241,13 +241,17 @@ const EMPTY_JAC = Matrix{Num}(undef, 0, 0) function invalidate_cache!(sys::AbstractSystem) if has_tgrad(sys) get_tgrad(sys)[] = EMPTY_TGRAD - elseif has_jac(sys) + end + if has_jac(sys) get_jac(sys)[] = EMPTY_JAC - elseif has_ctrl_jac(sys) + end + if has_ctrl_jac(sys) get_ctrl_jac(sys)[] = EMPTY_JAC - elseif has_Wfact(sys) + end + if has_Wfact(sys) get_Wfact(sys)[] = EMPTY_JAC - elseif has_Wfact_t(sys) + end + if has_Wfact_t(sys) get_Wfact_t(sys)[] = EMPTY_JAC end return sys From 9e7fa1468cfbc1c218c01f409ad91e80cba2e4dc Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Wed, 13 Apr 2022 01:31:19 -0400 Subject: [PATCH 15/19] Why not start all over again? --- src/systems/abstractsystem.jl | 16 +- src/systems/connectors.jl | 304 +++++++++++++++++++--------------- 2 files changed, 177 insertions(+), 143 deletions(-) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 0dda0df4a2..b3ba561469 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -418,9 +418,9 @@ function namespace_equations(sys::AbstractSystem) map(eq->namespace_equation(eq, sys), eqs) end -function namespace_equation(eq::Equation, sys) - _lhs = namespace_expr(eq.lhs, sys) - _rhs = namespace_expr(eq.rhs, sys) +function namespace_equation(eq::Equation, sys, n=nameof(sys)) + _lhs = namespace_expr(eq.lhs, sys, n) + _rhs = namespace_expr(eq.rhs, sys, n) _lhs ~ _rhs end @@ -430,22 +430,22 @@ function namespace_assignment(eq::Assignment, sys) Assignment(_lhs, _rhs) end -function namespace_expr(O, sys) where {T} +function namespace_expr(O, sys, n=nameof(sys)) where {T} ivs = independent_variables(sys) O = unwrap(O) if any(isequal(O), ivs) return O elseif isvariable(O) - renamespace(sys, O) + renamespace(n, O) elseif istree(O) - renamed = map(a->namespace_expr(a, sys), arguments(O)) + renamed = map(a->namespace_expr(a, sys, n), arguments(O)) if symtype(operation(O)) <: FnType - renamespace(sys, O) + renamespace(n, O) else similarterm(O, operation(O), renamed) end elseif O isa Array - map(Base.Fix2(namespace_expr, sys), O) + map(o->namespace_expr(o, sys, n), O) else O end diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index fa1b94764d..2b66077d7b 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -285,26 +285,27 @@ function connection2set!(connectionsets, namespace, ss, isouter) end end -generate_connection_set(sys::AbstractSystem) = (connectionsets = ConnectionSet[]; (generate_connection_set!(connectionsets, sys::AbstractSystem), merge(connectionsets))) +function generate_connection_set(sys::AbstractSystem) + connectionsets = ConnectionSet[] + sys = generate_connection_set!(connectionsets, sys) + sys, merge(connectionsets) +end + function generate_connection_set!(connectionsets, sys::AbstractSystem, namespace=nothing) + @show namespace + subsys = get_systems(sys) - # no connectors if there are no subsystems - isempty(subsys) && return sys isouter = generate_isouter(sys) eqs′ = get_eqs(sys) - #eqs = Equation[] + eqs = Equation[] - #instream_eqs = Equation[] - #instream_exprs = [] cts = [] # connections for eq in eqs′ if eq.lhs isa Connection push!(cts, get_systems(eq.rhs)) - #elseif collect_instream!(instream_exprs, eq) - # push!(instream_eqs, eq) else - #push!(eqs, eq) # split connections and equations + push!(eqs, eq) # split connections and equations end end @@ -321,16 +322,13 @@ function generate_connection_set!(connectionsets, sys::AbstractSystem, namespace end end - # if there are no connections, we are done - isempty(cts) && return sys - for ct in cts connection2set!(connectionsets, namespace, ct, isouter) end # pre order traversal @set! sys.systems = map(s->generate_connection_set!(connectionsets, s, renamespace(namespace, nameof(s))), subsys) - #@set! sys.eqs = eqs + @set! sys.eqs = eqs end function Base.merge(csets::AbstractVector{<:ConnectionSet}) @@ -380,105 +378,120 @@ end function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10) sys, csets = generate_connection_set(sys) - ceqs, _ = generate_connection_equations_and_stream_connections(csets) - sys = _expand_connections(sys; debug=debug, tol=tol) - @set! sys.eqs = [equations(sys); ceqs] + ceqs, instream_csets = generate_connection_equations_and_stream_connections(csets) + additional_eqs = Equation[] + _sys = expand_instream2(instream_csets, sys; debug=debug, tol=tol) + Main._a[] = _sys + sys = flatten(sys) + @set! sys.eqs = [equations(_sys); ceqs; additional_eqs] end -function _expand_connections(sys::AbstractSystem; debug=false, tol=1e-10, - rename=Ref{Union{Nothing,Tuple{Symbol,Int}}}(nothing), stream_connects=[]) - subsys = get_systems(sys) - isempty(subsys) && return sys +function unnamespace(root, namespace) + root === nothing && return namespace + root = string(root) + namespace = string(namespace) + if length(namespace) > length(root) + @assert root == namespace[1:length(root)] + Symbol(namespace[nextind(namespace, length(root)):end]) + else + @assert root == namespace + nothing + end +end +function expand_instream2(csets::AbstractVector{<:ConnectionSet}, sys::AbstractSystem, namespace=nothing, prevnamespace=nothing; debug=false, tol=1e-8) + subsys = get_systems(sys) + # no connectors if there are no subsystems + #isempty(subsys) && return sys # post order traversal - @set! sys.systems = map(s->_expand_connections(s, debug=debug, tol=tol, - rename=rename, stream_connects=stream_connects), subsys) + @set! sys.systems = map(s->expand_instream2(csets, s, renamespace(namespace, nameof(s)), namespace; debug, tol), subsys) + subsys = get_systems(sys) + @info "Expanding" namespace - isouter = generate_isouter(sys) - sys = flatten(sys) - eqs′ = get_eqs(sys) + sub = Dict() eqs = Equation[] instream_eqs = Equation[] - instream_exprs = [] - cts = [] # connections - for eq in eqs′ - if eq.lhs isa Connection - push!(cts, get_systems(eq.rhs)) - elseif collect_instream!(instream_exprs, eq) - push!(instream_eqs, eq) - else - push!(eqs, eq) # split connections and equations + instream_exprs = Set() + for s in subsys + seqs = map(Base.Fix2(namespace_equation, s), get_eqs(s)) + for eq in seqs + if collect_instream!(instream_exprs, eq) + push!(instream_eqs, eq) + else + push!(eqs, eq) + end end + end - # if there are no connections, we are done - isempty(cts) && return sys - - sys2idx = Dict{Symbol,Int}() # system (name) to n-th connect statement - narg_connects = Connection[] - for (i, syss) in enumerate(cts) - # find intersecting connections - exclude = 0 # exclude the intersecting system - idx = 0 # idx of narg_connects - for (j, s) in enumerate(syss) - idx′ = get(sys2idx, nameof(s), nothing) - idx′ === nothing && continue - idx = idx′ - exclude = j - end - if exclude == 0 - outers = [] - inners = [] - for s in syss - isouter(s) ? push!(outers, s) : push!(inners, s) + for ex in instream_exprs + cset, idx_in_set, sv = get_cset_sv(namespace, ex, csets) + + n_inners = n_outers = 0 + for (i, e) in enumerate(cset) + if e.isouter + n_outers += 1 + else + n_inners += 1 end - push!(narg_connects, Connection(outers=outers, inners=inners)) - for s in syss - sys2idx[nameof(s)] = length(narg_connects) + end + @show ex idx_in_set ConnectionSet(cset) + @show n_inners, n_outers + if n_inners == 1 && n_outers == 0 + sub[ex] = sv + elseif n_inners == 2 && n_outers == 0 + other = idx_in_set == 1 ? 2 : 1 + sub[ex] = states(renamespace(unnamespace(namespace, cset[other].sys.namespace), cset[other].sys.sys), sv) + elseif n_inners == 1 && n_outers == 1 + if !cset[idx_in_set].isouter + other = idx_in_set == 1 ? 2 : 1 + outerstream = states(renamespace(unnamespace(namespace, cset[other].sys.namespace), cset[other].sys.sys), sv) + sub[ex] = instream(outerstream) end else - # fuse intersecting connections - for (j, s) in enumerate(syss); j == exclude && continue - sys2idx[nameof(s)] = idx - c = narg_connects[idx] - isouter(s) ? push!(c.outers, s) : push!(c.inners, s) + if !cset[idx_in_set].isouter + fv = flowvar(first(connectors)) + # mj.c.m_flow + innerfvs = [unwrap(states(s, fv)) for (j, s) in enumerate(cset) if j != idx_in_set && !s.isouter] + innersvs = [unwrap(states(s, sv)) for (j, s) in enumerate(cset) if j != idx_in_set && !s.isouter] + # ck.m_flow + outerfvs = [unwrap(states(s, fv)) for s in cset if s.isouter] + outersvs = [unwrap(states(s, sv)) for s in cset if s.isouter] + + sub[ex_n] = term(instream_rt, Val(length(innerfvs)), Val(length(outerfvs)), innerfvs..., innersvs..., outerfvs..., outersvs...) end end end - isempty(narg_connects) && error("Unreachable reached. Please file an issue.") - - @set! sys.connections = narg_connects - - # Bad things happen when there are more than one intersections - for c in narg_connects - @unpack outers, inners = c - len = length(outers) + length(inners) - allconnectors = Iterators.flatten((outers, inners)) - dups = find_duplicates(nameof(c) for c in allconnectors) - length(dups) == 0 || error("Connection([$(join(map(nameof, allconnectors), ", "))]) has duplicated connections: [$(join(collect(dups), ", "))].") - end + display(instream_exprs) + display(sub) + @set! sys.systems = [] + @set! sys.eqs = [get_eqs(sys); eqs; substitute(instream_eqs, sub)] + sys +end - # stream variables - if rename[] !== nothing - name, depth = rename[] - nc = length(stream_connects) - for i in nc-depth+1:nc - stream_connects[i] = renamespace(:room, stream_connects[i]) +function get_cset_sv(namespace, ex, csets) + ns_sv = only(arguments(ex)) + full_name_sv = renamespace(namespace, ns_sv) + + cidx = -1 + idx_in_set = -1 + sv = ns_sv + for (i, c) in enumerate(csets) + for (j, v) in enumerate(c.set) + if isequal(namespaced_var(v), full_name_sv) + cidx = i + idx_in_set = j + sv = v.v + end end end - nsc = 0 - for c in narg_connects - if isstreamconnection(c) - push!(stream_connects, c) - nsc += 1 - end + cidx < 0 && error("$ns_sv is not a variable inside stream connectors") + cset = csets[cidx].set + if namespace != first(cset).sys.namespace + cset = map(c->@set(c.isouter = false), cset) end - rename[] = nameof(sys), nsc - instream_eqs, additional_eqs = expand_instream(instream_eqs, instream_exprs, stream_connects; debug=debug, tol=tol) - - @set! sys.eqs = [eqs; instream_eqs; additional_eqs] - return sys + cset, idx_in_set, sv end @generated function _instream_split(::Val{inner_n}, ::Val{outer_n}, vars::NTuple{N,Any}) where {inner_n, outer_n, N} @@ -542,65 +555,85 @@ function instream_rt(ins::Val{inner_n}, outs::Val{outer_n}, vars::Vararg{Any,N}) end SymbolicUtils.promote_symtype(::typeof(instream_rt), ::Vararg) = Real -function expand_instream2(csets::AbstractVector{<:ConnectionSet}, sys::AbstractSystem, namespace=nothing, debug=false, tol=1e-8) +#= +function expand_instream2!(additional_eqs, csets::AbstractVector{<:ConnectionSet}, ins, sys::AbstractSystem, namespace=nothing, prevnamespace=nothing; debug=false, tol=1e-8) subsys = get_systems(sys) # no connectors if there are no subsystems - isempty(subsys) && return sys + isempty(subsys) && return # post order traversal - @set! sys.systems = map(s->expand_instream2(csets, s, renamespace(namespace, nameof(s)), debug, tol), subsys) - - eqs′ = get_eqs(sys) - eqs = Equation[] - instream_eqs = Equation[] - instream_exprs = [] - for eq in eqs′ - if collect_instream!(instream_exprs, eq) - push!(instream_eqs, eq) - else - push!(eqs, eq) # split connections and equations - end + for s in subsys + expand_instream2!(additional_eqs, csets, ins, s, renamespace(namespace, nameof(s)), namespace; debug, tol) end + instream_eqs, instream_exprs = ins[namespace] + sys = flatten(sys) sub = Dict() - for ex in instream_exprs - sv = only(arguments(ex)) - full_name_sv = renamespace(namespace, sv) + dels = Int[] + for (k, ex) in enumerate(instream_exprs) + ex_n = namespace_expr(ex, sys, namespace) + ns_sv = only(arguments(ex)) + full_name_sv = renamespace(namespace, ns_sv) + @show full_name_sv - #cidx = findfirst(c->any(v->, ), ) cidx = -1 idx_in_set = -1 + sv = ns_sv for (i, c) in enumerate(csets) for (j, v) in enumerate(c.set) if isequal(namespaced_var(v), full_name_sv) cidx = i idx_in_set = j + sv = v.v end end end - cidx < 0 && error("$sv is not a variable inside stream connectors") - cset = csets[cidx].set - - connectors = Vector{Any}(undef, length(cset)) - n_inners = n_outers = 0 - for (i, e) in enumerate(cset) - connectors[i] = e.sys.sys - if e.isouter - n_outers += 1 - else - n_inners += 1 + #cidx < 0 && error("$ns_sv is not a variable inside stream connectors") + if cidx > 0 + cset = csets[cidx].set + + connectors = Vector{Any}(undef, length(cset)) + n_inners = n_outers = 0 + for (i, e) in enumerate(cset) + connectors[i] = e.sys.sys + if e.isouter + n_outers += 1 + else + n_inners += 1 + end + end + #@assert all(s->first(cset).sys.namespace == s.sys.namespace, cset) + if first(cset).sys.namespace != prevnamespace + cset = map(c->@set(c.isouter = false), cset) end + @show prevnamespace + @show n_inners, n_outers + @show nameof(sys) ConnectionSet(cset) + else + n_inners = 1 + n_outers = 0 end if n_inners == 1 && n_outers == 0 - sub[ex] = sv + @show namespace, ex, sv + sub[ex_n] = renamespace(renamespace(namespace, nameof(sys)), sv) elseif n_inners == 2 && n_outers == 0 other = idx_in_set == 1 ? 2 : 1 - sub[ex] = states(cset[other], sv) + sub[ex_n] = states(cset[other], sv) elseif n_inners == 1 && n_outers == 1 if !cset[idx_in_set].isouter other = idx_in_set == 1 ? 2 : 1 - outerstream = states(cset[other], sv) - sub[ex] = outerstream + outerstream = instream(states(cset[other].sys.sys, sv)) + @show outerstream + ns = nameof(sys) + ex = namespace_expr(ex, sys, ns) + # TODO: write a mapping from eqs to exprs + push!(dels, k) + eq = substitute(namespace_equation(instream_eqs[k], sys, ns), Dict(ex=>outerstream)) + previnstream_eqs, previnstream_exprs = ins[prevnamespace] + push!(previnstream_eqs, eq) + push!(previnstream_exprs, outerstream) + #outerstream = states(cset[other], sv) + #substitute() end else if !cset[idx_in_set].isouter @@ -612,13 +645,14 @@ function expand_instream2(csets::AbstractVector{<:ConnectionSet}, sys::AbstractS outerfvs = [unwrap(states(s, fv)) for s in cset if s.isouter] outersvs = [unwrap(states(s, sv)) for s in cset if s.isouter] - sub[ex] = term(instream_rt, Val(length(innerfvs)), Val(length(outerfvs)), innerfvs..., innersvs..., outerfvs..., outersvs...) + sub[ex_n] = term(instream_rt, Val(length(innerfvs)), Val(length(outerfvs)), innerfvs..., innersvs..., outerfvs..., outersvs...) end end end + instream_eqs = deleteat!(copy(instream_eqs), dels) + instream_eqs = map(eq->substitute(namespace_equation(eq, sys, namespace), sub), instream_eqs) # additional equations - additional_eqs = Equation[] csets = filter(cset->any(e->e.sys.namespace === namespace, cset.set), csets) for cset′ in csets cset = cset′.set @@ -676,7 +710,6 @@ function expand_instream2(csets::AbstractVector{<:ConnectionSet}, sys::AbstractS end end - instream_eqs = map(Base.Fix2(substitute, sub), instream_eqs) if debug println("===========BEGIN=============") println("Expanded equations:") @@ -691,9 +724,10 @@ function expand_instream2(csets::AbstractVector{<:ConnectionSet}, sys::AbstractS end println("============END==============") end - - @set! sys.eqs = [eqs; instream_eqs; additional_eqs] + append!(additional_eqs, instream_eqs) + return end +=# function expand_instream(instream_eqs, instream_exprs, connects; debug=false, tol) sub = Dict() @@ -727,19 +761,19 @@ function expand_instream(instream_eqs, instream_exprs, connects; debug=false, to # recommended: if n_inners == 1 && n_outers == 0 connector_name === only(inner_names) || error("$var is not in any stream connector of $(nameof(ogsys))") - sub[ex] = var + sub[ex_n] = var elseif n_inners == 2 && n_outers == 0 connector_name in inner_names || error("$var is not in any stream connector of $(nameof(ogsys))") idx = findfirst(c->nameof(c) === connector_name, inner_sc) other = idx == 1 ? 2 : 1 - sub[ex] = states(inner_sc[other], sv) + sub[ex_n] = states(inner_sc[other], sv) elseif n_inners == 1 && n_outers == 1 isinner = connector_name === only(inner_names) isouter = connector_name === only(outer_names) (isinner || isouter) || error("$var is not in any stream connector of $(nameof(ogsys))") if isinner outerstream = states(only(outer_sc), sv) # c_1.h_outflow - sub[ex] = outerstream + sub[ex_n] = outerstream end else fv = flowvar(first(connectors)) @@ -752,7 +786,7 @@ function expand_instream(instream_eqs, instream_exprs, connects; debug=false, to outerfvs = [unwrap(states(s, fv)) for s in outer_sc] outersvs = [instream(states(s, sv)) for s in outer_sc] - sub[ex] = term(instream_rt, Val(length(innerfvs)), Val(length(outerfvs)), innerfvs..., innersvs..., outerfvs..., outersvs...) + sub[ex_n] = term(instream_rt, Val(length(innerfvs)), Val(length(outerfvs)), innerfvs..., innersvs..., outerfvs..., outersvs...) #= si = isempty(outer_sc) ? 0 : sum(s->max(states(s, fv), 0), outer_sc) for j in 1:n_inners; j == i && continue @@ -774,7 +808,7 @@ function expand_instream(instream_eqs, instream_exprs, connects; debug=false, to den += tmp num += tmp * instream(states(outer_sc[k], sv)) end - sub[ex] = mydiv(num, den) + sub[ex_n] = mydiv(num, den) =# end end @@ -847,5 +881,5 @@ function expand_instream(instream_eqs, instream_exprs, connects; debug=false, to end println("============END==============") end - return instream_eqs, additional_eqs + return end From 039b86c334db64f598810b0789d806d70258a0be Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Wed, 13 Apr 2022 12:48:56 -0400 Subject: [PATCH 16/19] Stream connector rewrite --- src/systems/connectors.jl | 83 ++++++++++++++++++++++++++++++++++----- 1 file changed, 74 insertions(+), 9 deletions(-) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index 2b66077d7b..986bdf8244 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -441,35 +441,100 @@ function expand_instream2(csets::AbstractVector{<:ConnectionSet}, sys::AbstractS sub[ex] = sv elseif n_inners == 2 && n_outers == 0 other = idx_in_set == 1 ? 2 : 1 - sub[ex] = states(renamespace(unnamespace(namespace, cset[other].sys.namespace), cset[other].sys.sys), sv) + sub[ex] = get_current_var(namespace, cset[other], sv) elseif n_inners == 1 && n_outers == 1 if !cset[idx_in_set].isouter other = idx_in_set == 1 ? 2 : 1 - outerstream = states(renamespace(unnamespace(namespace, cset[other].sys.namespace), cset[other].sys.sys), sv) + outerstream = get_current_var(namespace, cset[other], sv) sub[ex] = instream(outerstream) end else if !cset[idx_in_set].isouter - fv = flowvar(first(connectors)) + fv = flowvar(first(cset).sys.sys) # mj.c.m_flow - innerfvs = [unwrap(states(s, fv)) for (j, s) in enumerate(cset) if j != idx_in_set && !s.isouter] - innersvs = [unwrap(states(s, sv)) for (j, s) in enumerate(cset) if j != idx_in_set && !s.isouter] + innerfvs = [get_current_var(namespace, s, fv) for (j, s) in enumerate(cset) if j != idx_in_set && !s.isouter] + innersvs = [get_current_var(namespace, s, sv) for (j, s) in enumerate(cset) if j != idx_in_set && !s.isouter] # ck.m_flow - outerfvs = [unwrap(states(s, fv)) for s in cset if s.isouter] - outersvs = [unwrap(states(s, sv)) for s in cset if s.isouter] + outerfvs = [get_current_var(namespace, s, fv) for s in cset if s.isouter] + outersvs = [get_current_var(namespace, s, sv) for s in cset if s.isouter] - sub[ex_n] = term(instream_rt, Val(length(innerfvs)), Val(length(outerfvs)), innerfvs..., innersvs..., outerfvs..., outersvs...) + sub[ex] = term(instream_rt, Val(length(innerfvs)), Val(length(outerfvs)), innerfvs..., innersvs..., outerfvs..., outersvs...) + end + end + end + + # additional equations + additional_eqs = Equation[] + csets = filter(cset->any(e->e.sys.namespace === namespace, cset.set), csets) + @show csets + for cset′ in csets + cset = cset′.set + connectors = Vector{Any}(undef, length(cset)) + n_inners = n_outers = 0 + for (i, e) in enumerate(cset) + connectors[i] = e.sys.sys + if e.isouter + n_outers += 1 + else + n_inners += 1 + end + end + iszero(n_outers) && continue + connector_representative = first(cset).sys.sys + fv = flowvar(connector_representative) + sv = first(cset).v + vtype = get_connection_type(sv) + vtype === Stream || continue + if n_inners == 1 && n_outers == 1 + push!(additional_eqs, @show states(cset[1].sys.sys, sv) ~ states(cset[2].sys.sys, sv)) + elseif n_inners == 0 && n_outers == 2 + # we don't expand `instream` in this case. + v1 = states(cset[1].sys.sys, sv) + v2 = states(cset[2].sys.sys, sv) + push!(additional_eqs, v1 ~ instream(v2)) + push!(additional_eqs, v2 ~ instream(v1)) + else + sq = 0 + s_inners = (s for s in cset if !s.isouter) + s_outers = (s for s in cset if s.isouter) + for (q, oscq) in enumerate(s_outers) + sq += sum(s->max(-states(s, fv), 0), s_inners) + for (k, s) in enumerate(s_outers); k == q && continue + f = states(s.sys.sys, fv) + sq += max(f, 0) + end + + num = 0 + den = 0 + for s in s_inners + f = states(s.sys.sys, fv) + tmp = positivemax(-f, sq; tol=tol) + den += tmp + num += tmp * states(s.sys.sys, sv) + end + for (k, s) in enumerate(s_outers); k == q && continue + f = states(s.sys.sys, fv) + tmp = positivemax(f, sq; tol=tol) + den += tmp + num += tmp * instream(states(s.sys.sys, sv)) + end + push!(additional_eqs, states(oscq.sys.sys, sv) ~ num / den) end end end + @show additional_eqs display(instream_exprs) display(sub) @set! sys.systems = [] - @set! sys.eqs = [get_eqs(sys); eqs; substitute(instream_eqs, sub)] + @set! sys.eqs = [get_eqs(sys); eqs; substitute(instream_eqs, sub); additional_eqs] sys end +function get_current_var(namespace, cele, sv) + states(renamespace(unnamespace(namespace, cele.sys.namespace), cele.sys.sys), sv) +end + function get_cset_sv(namespace, ex, csets) ns_sv = only(arguments(ex)) full_name_sv = renamespace(namespace, ns_sv) From 9c166ccdea480598a2ffcafb63841f7dd353550e Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Wed, 13 Apr 2022 13:11:54 -0400 Subject: [PATCH 17/19] Clean up --- src/systems/connectors.jl | 405 ++------------------------------- src/systems/systemstructure.jl | 8 +- 2 files changed, 22 insertions(+), 391 deletions(-) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index 986bdf8244..b01216fe4f 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -93,54 +93,6 @@ function connect(sys1::AbstractSystem, sys2::AbstractSystem, syss::AbstractSyste Equation(Connection(), Connection(syss)) # the RHS are connected systems end -# the actual `connect`. -function connect(c::Connection; check=true) - @unpack inners, outers = c - - cnts = Iterators.flatten((inners, outers)) - fs, ss = Iterators.peel(cnts) - splitting_idx = length(inners) # anything after the splitting_idx is outer. - first_sts = get_states(fs) - first_sts_set = Set(getname.(first_sts)) - for sys in ss - current_sts = getname.(get_states(sys)) - Set(current_sts) == first_sts_set || error("$(nameof(sys)) ($current_sts) doesn't match the connection type of $(nameof(fs)) ($first_sts).") - end - - seen = Set() - ceqs = Equation[] - for s in first_sts - name = getname(s) - fix_val = getproperty(fs, name) # representative - fix_val in seen && continue - push!(seen, fix_val) - - vtype = get_connection_type(fix_val) - vtype === Stream && continue - - if vtype === Flow - for j in eachindex(fix_val) - rhs = 0 - for (i, c) in enumerate(cnts) - isinner = i <= splitting_idx - var = getproperty(c, name) - rhs += isinner ? var[j] : -var[j] - end - push!(ceqs, 0 ~ rhs) - end - else # Equality - for c in ss - var = getproperty(c, name) - for (i, v) in enumerate(var) - push!(ceqs, fix_val[i] ~ v) - end - end - end - end - - return ceqs -end - instream(a) = term(instream, unwrap(a), type=symtype(a)) SymbolicUtils.promote_symtype(::typeof(instream), _) = Real @@ -292,8 +244,6 @@ function generate_connection_set(sys::AbstractSystem) end function generate_connection_set!(connectionsets, sys::AbstractSystem, namespace=nothing) - @show namespace - subsys = get_systems(sys) isouter = generate_isouter(sys) @@ -381,7 +331,6 @@ function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10) ceqs, instream_csets = generate_connection_equations_and_stream_connections(csets) additional_eqs = Equation[] _sys = expand_instream2(instream_csets, sys; debug=debug, tol=tol) - Main._a[] = _sys sys = flatten(sys) @set! sys.eqs = [equations(_sys); ceqs; additional_eqs] end @@ -406,7 +355,10 @@ function expand_instream2(csets::AbstractVector{<:ConnectionSet}, sys::AbstractS # post order traversal @set! sys.systems = map(s->expand_instream2(csets, s, renamespace(namespace, nameof(s)), namespace; debug, tol), subsys) subsys = get_systems(sys) - @info "Expanding" namespace + + if debug + @info "Expanding" namespace + end sub = Dict() eqs = Equation[] @@ -435,8 +387,10 @@ function expand_instream2(csets::AbstractVector{<:ConnectionSet}, sys::AbstractS n_inners += 1 end end - @show ex idx_in_set ConnectionSet(cset) - @show n_inners, n_outers + if debug + @show ex idx_in_set ConnectionSet(cset) + @show n_inners, n_outers + end if n_inners == 1 && n_outers == 0 sub[ex] = sv elseif n_inners == 2 && n_outers == 0 @@ -466,7 +420,6 @@ function expand_instream2(csets::AbstractVector{<:ConnectionSet}, sys::AbstractS # additional equations additional_eqs = Equation[] csets = filter(cset->any(e->e.sys.namespace === namespace, cset.set), csets) - @show csets for cset′ in csets cset = cset′.set connectors = Vector{Any}(undef, length(cset)) @@ -486,7 +439,7 @@ function expand_instream2(csets::AbstractVector{<:ConnectionSet}, sys::AbstractS vtype = get_connection_type(sv) vtype === Stream || continue if n_inners == 1 && n_outers == 1 - push!(additional_eqs, @show states(cset[1].sys.sys, sv) ~ states(cset[2].sys.sys, sv)) + push!(additional_eqs, states(cset[1].sys.sys, sv) ~ states(cset[2].sys.sys, sv)) elseif n_inners == 0 && n_outers == 2 # we don't expand `instream` in this case. v1 = states(cset[1].sys.sys, sv) @@ -522,10 +475,14 @@ function expand_instream2(csets::AbstractVector{<:ConnectionSet}, sys::AbstractS end end end - @show additional_eqs - display(instream_exprs) - display(sub) + if debug + @info "Additional equations" csets + @show additional_eqs + println("Substitutions") + display(sub) + end + @set! sys.systems = [] @set! sys.eqs = [get_eqs(sys); eqs; substitute(instream_eqs, sub); additional_eqs] sys @@ -559,6 +516,7 @@ function get_cset_sv(namespace, ex, csets) cset, idx_in_set, sv end +# instream runtime @generated function _instream_split(::Val{inner_n}, ::Val{outer_n}, vars::NTuple{N,Any}) where {inner_n, outer_n, N} #instream_rt(innerfvs..., innersvs..., outerfvs..., outersvs...) ret = Expr(:tuple) @@ -619,332 +577,3 @@ function instream_rt(ins::Val{inner_n}, outs::Val{outer_n}, vars::Vararg{Any,N}) =# end SymbolicUtils.promote_symtype(::typeof(instream_rt), ::Vararg) = Real - -#= -function expand_instream2!(additional_eqs, csets::AbstractVector{<:ConnectionSet}, ins, sys::AbstractSystem, namespace=nothing, prevnamespace=nothing; debug=false, tol=1e-8) - subsys = get_systems(sys) - # no connectors if there are no subsystems - isempty(subsys) && return - # post order traversal - for s in subsys - expand_instream2!(additional_eqs, csets, ins, s, renamespace(namespace, nameof(s)), namespace; debug, tol) - end - - instream_eqs, instream_exprs = ins[namespace] - sys = flatten(sys) - - sub = Dict() - dels = Int[] - for (k, ex) in enumerate(instream_exprs) - ex_n = namespace_expr(ex, sys, namespace) - ns_sv = only(arguments(ex)) - full_name_sv = renamespace(namespace, ns_sv) - @show full_name_sv - - cidx = -1 - idx_in_set = -1 - sv = ns_sv - for (i, c) in enumerate(csets) - for (j, v) in enumerate(c.set) - if isequal(namespaced_var(v), full_name_sv) - cidx = i - idx_in_set = j - sv = v.v - end - end - end - #cidx < 0 && error("$ns_sv is not a variable inside stream connectors") - if cidx > 0 - cset = csets[cidx].set - - connectors = Vector{Any}(undef, length(cset)) - n_inners = n_outers = 0 - for (i, e) in enumerate(cset) - connectors[i] = e.sys.sys - if e.isouter - n_outers += 1 - else - n_inners += 1 - end - end - #@assert all(s->first(cset).sys.namespace == s.sys.namespace, cset) - if first(cset).sys.namespace != prevnamespace - cset = map(c->@set(c.isouter = false), cset) - end - @show prevnamespace - @show n_inners, n_outers - @show nameof(sys) ConnectionSet(cset) - else - n_inners = 1 - n_outers = 0 - end - if n_inners == 1 && n_outers == 0 - @show namespace, ex, sv - sub[ex_n] = renamespace(renamespace(namespace, nameof(sys)), sv) - elseif n_inners == 2 && n_outers == 0 - other = idx_in_set == 1 ? 2 : 1 - sub[ex_n] = states(cset[other], sv) - elseif n_inners == 1 && n_outers == 1 - if !cset[idx_in_set].isouter - other = idx_in_set == 1 ? 2 : 1 - outerstream = instream(states(cset[other].sys.sys, sv)) - @show outerstream - ns = nameof(sys) - ex = namespace_expr(ex, sys, ns) - # TODO: write a mapping from eqs to exprs - push!(dels, k) - eq = substitute(namespace_equation(instream_eqs[k], sys, ns), Dict(ex=>outerstream)) - previnstream_eqs, previnstream_exprs = ins[prevnamespace] - push!(previnstream_eqs, eq) - push!(previnstream_exprs, outerstream) - #outerstream = states(cset[other], sv) - #substitute() - end - else - if !cset[idx_in_set].isouter - fv = flowvar(first(connectors)) - # mj.c.m_flow - innerfvs = [unwrap(states(s, fv)) for (j, s) in enumerate(cset) if j != idx_in_set && !s.isouter] - innersvs = [unwrap(states(s, sv)) for (j, s) in enumerate(cset) if j != idx_in_set && !s.isouter] - # ck.m_flow - outerfvs = [unwrap(states(s, fv)) for s in cset if s.isouter] - outersvs = [unwrap(states(s, sv)) for s in cset if s.isouter] - - sub[ex_n] = term(instream_rt, Val(length(innerfvs)), Val(length(outerfvs)), innerfvs..., innersvs..., outerfvs..., outersvs...) - end - end - end - instream_eqs = deleteat!(copy(instream_eqs), dels) - instream_eqs = map(eq->substitute(namespace_equation(eq, sys, namespace), sub), instream_eqs) - - # additional equations - csets = filter(cset->any(e->e.sys.namespace === namespace, cset.set), csets) - for cset′ in csets - cset = cset′.set - connectors = Vector{Any}(undef, length(cset)) - n_inners = n_outers = 0 - for (i, e) in enumerate(cset) - connectors[i] = e.sys.sys - if e.isouter - n_outers += 1 - else - n_inners += 1 - end - end - iszero(n_outers) && continue - connector_representative = first(cset).sys.sys - fv = flowvar(connector_representative) - sv = first(cset).v - vtype = get_connection_type(sv) - vtype === Stream || continue - if n_inners == 1 && n_outers == 1 - push!(additional_eqs, states(cset[1], sv) ~ states(cset[2], sv)) - elseif n_inners == 0 && n_outers == 2 - # we don't expand `instream` in this case. - v1 = states(cset[1], sv) - v2 = states(cset[2], sv) - push!(additional_eqs, v1 ~ instream(v2)) - push!(additional_eqs, v2 ~ instream(v1)) - else - sq = 0 - s_inners = (s for s in cset if !s.isouter) - s_outers = (s for s in cset if s.isouter) - for (q, oscq) in enumerate(s_outers) - sq += sum(s->max(-states(s, fv), 0), s_inners) - for (k, s) in enumerate(s_outers); k == q && continue - f = states(s, fv) - sq += max(f, 0) - end - - num = 0 - den = 0 - for s in s_inners - f = states(s, fv) - tmp = positivemax(-f, sq; tol=tol) - den += tmp - num += tmp * states(s, sv) - end - for (k, s) in enumerate(s_outers); k == q && continue - f = states(s, fv) - tmp = positivemax(f, sq; tol=tol) - den += tmp - num += tmp * instream(states(s, sv)) - end - push!(additional_eqs, states(oscq, sv) ~ num / den) - end - end - end - - if debug - println("===========BEGIN=============") - println("Expanded equations:") - for eq in instream_eqs - print_with_indent(4, eq) - end - if !isempty(additional_eqs) - println("Additional equations:") - for eq in additional_eqs - print_with_indent(4, eq) - end - end - println("============END==============") - end - append!(additional_eqs, instream_eqs) - return -end -=# - -function expand_instream(instream_eqs, instream_exprs, connects; debug=false, tol) - sub = Dict() - for ex in instream_exprs - var = only(arguments(ex)) - connector_name, streamvar_name = split_sys_var(var) - - # find the connect - cidx = findfirst(c->connector_name in c, connects) - cidx === nothing && error("$var is not a variable inside stream connectors") - connect = connects[cidx] - connectors = Iterators.flatten((connect.inners, connect.outers)) - # stream variable - sv = getproperty(first(connectors), streamvar_name; namespace=false) - inner_sc = something(connect.inners, EMPTY_VEC) - outer_sc = something(connect.outers, EMPTY_VEC) - - n_outers = length(outer_sc) - n_inners = length(inner_sc) - outer_names = (nameof(s) for s in outer_sc) - inner_names = (nameof(s) for s in inner_sc) - if debug - println("Expanding: $ex") - isempty(inner_names) || println("Inner connectors: $(collect(inner_names))") - isempty(outer_names) || println("Outer connectors: $(collect(outer_names))") - end - - # expand `instream`s - # https://specification.modelica.org/v3.4/Ch15.html - # Based on the above requirements, the following implementation is - # recommended: - if n_inners == 1 && n_outers == 0 - connector_name === only(inner_names) || error("$var is not in any stream connector of $(nameof(ogsys))") - sub[ex_n] = var - elseif n_inners == 2 && n_outers == 0 - connector_name in inner_names || error("$var is not in any stream connector of $(nameof(ogsys))") - idx = findfirst(c->nameof(c) === connector_name, inner_sc) - other = idx == 1 ? 2 : 1 - sub[ex_n] = states(inner_sc[other], sv) - elseif n_inners == 1 && n_outers == 1 - isinner = connector_name === only(inner_names) - isouter = connector_name === only(outer_names) - (isinner || isouter) || error("$var is not in any stream connector of $(nameof(ogsys))") - if isinner - outerstream = states(only(outer_sc), sv) # c_1.h_outflow - sub[ex_n] = outerstream - end - else - fv = flowvar(first(connectors)) - i = findfirst(c->nameof(c) === connector_name, inner_sc) - if i !== nothing - # mj.c.m_flow - innerfvs = [unwrap(states(s, fv)) for (j, s) in enumerate(inner_sc) if j != i] - innersvs = [unwrap(states(s, sv)) for (j, s) in enumerate(inner_sc) if j != i] - # ck.m_flow - outerfvs = [unwrap(states(s, fv)) for s in outer_sc] - outersvs = [instream(states(s, sv)) for s in outer_sc] - - sub[ex_n] = term(instream_rt, Val(length(innerfvs)), Val(length(outerfvs)), innerfvs..., innersvs..., outerfvs..., outersvs...) - #= - si = isempty(outer_sc) ? 0 : sum(s->max(states(s, fv), 0), outer_sc) - for j in 1:n_inners; j == i && continue - f = states(inner_sc[j], fv) - si += max(-f, 0) - end - - num = 0 - den = 0 - for j in 1:n_inners; j == i && continue - f = states(inner_sc[j], fv) - tmp = positivemax(-f, si; tol=tol) - den += tmp - num += tmp * states(inner_sc[j], sv) - end - for k in 1:n_outers - f = states(outer_sc[k], fv) - tmp = positivemax(f, si; tol=tol) - den += tmp - num += tmp * instream(states(outer_sc[k], sv)) - end - sub[ex_n] = mydiv(num, den) - =# - end - end - end - - # additional equations - additional_eqs = Equation[] - for c in connects - outer_sc = something(c.outers, EMPTY_VEC) - isempty(outer_sc) && continue - inner_sc = c.inners - n_outers = length(outer_sc) - n_inners = length(inner_sc) - connector_representative = first(outer_sc) - fv = flowvar(connector_representative) - for sv in get_states(connector_representative) - vtype = get_connection_type(sv) - vtype === Stream || continue - if n_inners == 1 && n_outers == 1 - innerstream = states(only(inner_sc), sv) - outerstream = states(only(outer_sc), sv) - push!(additional_eqs, outerstream ~ innerstream) - elseif n_inners == 0 && n_outers == 2 - # we don't expand `instream` in this case. - v1 = states(outer_sc[1], sv) - v2 = states(outer_sc[2], sv) - push!(additional_eqs, v1 ~ instream(v2)) - push!(additional_eqs, v2 ~ instream(v1)) - else - sq = 0 - for q in 1:n_outers - sq += sum(s->max(-states(s, fv), 0), inner_sc) - for k in 1:n_outers; k == q && continue - f = states(outer_sc[k], fv) - sq += max(f, 0) - end - - num = 0 - den = 0 - for j in 1:n_inners - f = states(inner_sc[j], fv) - tmp = positivemax(-f, sq; tol=tol) - den += tmp - num += tmp * states(inner_sc[j], sv) - end - for k in 1:n_outers; k == q && continue - f = states(outer_sc[k], fv) - tmp = positivemax(f, sq; tol=tol) - den += tmp - num += tmp * instream(states(outer_sc[k], sv)) - end - push!(additional_eqs, states(outer_sc[q], sv) ~ num / den) - end - end - end - end - - instream_eqs = map(Base.Fix2(substitute, sub), instream_eqs) - if debug - println("===========BEGIN=============") - println("Expanded equations:") - for eq in instream_eqs - print_with_indent(4, eq) - end - if !isempty(additional_eqs) - println("Additional equations:") - for eq in additional_eqs - print_with_indent(4, eq) - end - end - println("============END==============") - end - return -end diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 51d6ac374c..50e8878457 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -299,10 +299,12 @@ function TearingState(sys; quick_cancel=false, check=true) # it could be that a variable appeared in the states, but never appeared # in the equations. algvaridx = get(var2idx, algvar, 0) - if algvaridx == 0 - check ? throw(InvalidSystemException("The system is missing an equation for $algvar.")) : return nothing + #if algvaridx == 0 + # check ? throw(InvalidSystemException("The system is missing an equation for $algvar.")) : return nothing + #end + if algvaridx != 0 + vartype[algvaridx] = ALGEBRAIC_VARIABLE end - vartype[algvaridx] = ALGEBRAIC_VARIABLE end graph = BipartiteGraph(neqs, nvars, Val(false)) From b7e4cefc02a4378c606d0323744396590f9b14b3 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Wed, 13 Apr 2022 14:57:04 -0400 Subject: [PATCH 18/19] More performance --- src/systems/connectors.jl | 39 +++++++++++++++++------- src/systems/diffeqs/odesystem.jl | 4 +-- src/systems/nonlinear/nonlinearsystem.jl | 4 +-- 3 files changed, 32 insertions(+), 15 deletions(-) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index b01216fe4f..bffa2af2b8 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -283,13 +283,32 @@ end function Base.merge(csets::AbstractVector{<:ConnectionSet}) mcsets = ConnectionSet[] - # FIXME: this is O(m n^3) + ele2idx = Dict{ConnectionElement,Int}() + cacheset = Set{ConnectionElement}() for cset in csets - idx = findfirst(mcset->any(s->any(z->z == s, cset.set), mcset.set), mcsets) + idx = nothing + for e in cset.set + idx = get(ele2idx, e, nothing) + idx !== nothing && break + end if idx === nothing push!(mcsets, cset) + for e in cset.set + ele2idx[e] = length(mcsets) + end else - union!(mcsets[idx].set, cset.set) + for e in mcsets[idx].set + push!(cacheset, e) + end + for e in cset.set + push!(cacheset, e) + end + empty!(mcsets[idx].set) + for e in cacheset + ele2idx[e] = idx + push!(mcsets[idx].set, e) + end + empty!(cacheset) end end mcsets @@ -330,8 +349,8 @@ function expand_connections(sys::AbstractSystem; debug=false, tol=1e-10) sys, csets = generate_connection_set(sys) ceqs, instream_csets = generate_connection_equations_and_stream_connections(csets) additional_eqs = Equation[] - _sys = expand_instream2(instream_csets, sys; debug=debug, tol=tol) - sys = flatten(sys) + _sys = expand_instream(instream_csets, sys; debug=debug, tol=tol) + sys = flatten(sys, true) @set! sys.eqs = [equations(_sys); ceqs; additional_eqs] end @@ -348,12 +367,10 @@ function unnamespace(root, namespace) end end -function expand_instream2(csets::AbstractVector{<:ConnectionSet}, sys::AbstractSystem, namespace=nothing, prevnamespace=nothing; debug=false, tol=1e-8) +function expand_instream(csets::AbstractVector{<:ConnectionSet}, sys::AbstractSystem, namespace=nothing, prevnamespace=nothing; debug=false, tol=1e-8) subsys = get_systems(sys) - # no connectors if there are no subsystems - #isempty(subsys) && return sys # post order traversal - @set! sys.systems = map(s->expand_instream2(csets, s, renamespace(namespace, nameof(s)), namespace; debug, tol), subsys) + @set! sys.systems = map(s->expand_instream(csets, s, renamespace(namespace, nameof(s)), namespace; debug, tol), subsys) subsys = get_systems(sys) if debug @@ -365,8 +382,8 @@ function expand_instream2(csets::AbstractVector{<:ConnectionSet}, sys::AbstractS instream_eqs = Equation[] instream_exprs = Set() for s in subsys - seqs = map(Base.Fix2(namespace_equation, s), get_eqs(s)) - for eq in seqs + for eq in get_eqs(s) + eq = namespace_equation(eq, s) if collect_instream!(instream_exprs, eq) push!(instream_eqs, eq) else diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index f845baad55..60707f2ed5 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -229,13 +229,13 @@ function Base.:(==)(sys1::ODESystem, sys2::ODESystem) all(s1 == s2 for (s1, s2) in zip(get_systems(sys1), get_systems(sys2))) end -function flatten(sys::ODESystem) +function flatten(sys::ODESystem, noeqs=false) systems = get_systems(sys) if isempty(systems) return sys else return ODESystem( - equations(sys), + noeqs ? Equation[] : equations(sys), get_iv(sys), states(sys), parameters(sys), diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index 38f3d4a497..01621d7125 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -335,13 +335,13 @@ function NonlinearProblemExpr{iip}(sys::NonlinearSystem,u0map, !linenumbers ? striplines(ex) : ex end -function flatten(sys::NonlinearSystem) +function flatten(sys::NonlinearSystem, noeqs=false) systems = get_systems(sys) if isempty(systems) return sys else return NonlinearSystem( - equations(sys), + noeqs ? Equation[] : equations(sys), states(sys), parameters(sys), observed=observed(sys), From 4dce14fa6800249a36c580f35c4867ac50137b40 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Wed, 13 Apr 2022 15:55:28 -0400 Subject: [PATCH 19/19] Fix test --- src/systems/abstractsystem.jl | 2 +- test/components.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index b3ba561469..d71b8189bb 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -509,7 +509,7 @@ for f in [:states, :parameters] @eval $f(sys::AbstractSystem, vs::AbstractArray) = map(v->$f(sys, v), vs) end -flatten(sys::AbstractSystem) = sys +flatten(sys::AbstractSystem, args...) = sys function equations(sys::ModelingToolkit.AbstractSystem) eqs = get_eqs(sys) diff --git a/test/components.jl b/test/components.jl index e0311aca5f..25ebacd16f 100644 --- a/test/components.jl +++ b/test/components.jl @@ -184,4 +184,4 @@ function Circuit(;name) end @named foo = Circuit() -@test_broken structural_simplify(foo) isa AbstractSystem +@test structural_simplify(foo) isa ModelingToolkit.AbstractSystem