Merge branch 'master' into myb/strong_alias

Author: YingboMa

src/bipartite_graph.jl

@@ -191,14 +191,49 @@ end
 # Matrix whose only purpose is to pretty-print the bipartite graph
 struct BipartiteAdjacencyList
     u::Union{Vector{Int}, Nothing}
+    highligh_u::Union{Set{Int}, Nothing}
+    match::Union{Int, Unassigned}
 end
+function BipartiteAdjacencyList(u::Union{Vector{Int}, Nothing})
+    BipartiteAdjacencyList(u, nothing, unassigned)
+end
+
+struct HighlightInt
+    i::Int
+    highlight::Union{Symbol, Nothing}
+end
+Base.typeinfo_implicit(::Type{HighlightInt}) = true
+
+function Base.show(io::IO, hi::HighlightInt)
+    if hi.highlight !== nothing
+        printstyled(io, hi.i, color = hi.highlight)
+    else
+        print(io, hi.i)
+    end
+end
+
 function Base.show(io::IO, l::BipartiteAdjacencyList)
     if l.u === nothing
         printstyled(io, '⋅', color = :light_black)
     elseif isempty(l.u)
         printstyled(io, '∅', color = :light_black)
-    else
+    elseif l.highligh_u === nothing
         print(io, l.u)
+    else
+        function choose_color(i)
+            i in l.highligh_u ? (i == l.match ? :light_yellow : :green) :
+            (i == l.match ? :yellow : nothing)
+        end
+        if !isempty(setdiff(l.highligh_u, l.u))
+            # Only for debugging, shouldn't happen in practice
+            print(io, map(union(l.u, l.highligh_u)) do i
+                      HighlightInt(i, !(i in l.u) ? :light_red : choose_color(i))
+                  end)
+        else
+            print(io, map(l.u) do i
+                      HighlightInt(i, choose_color(i))
+                  end)
+        end
     end
 end
 

src/inputoutput.jl

@@ -218,7 +218,8 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
     f, dvs, ps
 end
 
-function inputs_to_parameters!(state::TransformationState, check_bound = true)
+function inputs_to_parameters!(state::TransformationState, io)
+    check_bound = io === nothing
     @unpack structure, fullvars, sys = state
     @unpack var_to_diff, graph, solvable_graph = structure
     @assert solvable_graph === nothing
@@ -274,6 +275,17 @@ function inputs_to_parameters!(state::TransformationState, check_bound = true)
     @set! sys.eqs = map(Base.Fix2(substitute, input_to_parameters), equations(sys))
     @set! sys.states = setdiff(states(sys), keys(input_to_parameters))
     ps = parameters(sys)
+
+    if io !== nothing
+        # Change order of new parameters to correspond to user-provided order in argument `inputs`
+        d = Dict{Any, Int}()
+        for (i, inp) in enumerate(new_parameters)
+            d[inp] = i
+        end
+        permutation = [d[i] for i in io.inputs]
+        new_parameters = new_parameters[permutation]
+    end
+
     @set! sys.ps = [ps; new_parameters]
 
     @set! state.sys = sys

src/systems/abstractsystem.jl

@@ -968,7 +968,7 @@ function structural_simplify(sys::AbstractSystem, io = nothing; simplify = false
     state = TearingState(sys)
     has_io = io !== nothing
     has_io && markio!(state, io...)
-    state, input_idxs = inputs_to_parameters!(state, !has_io)
+    state, input_idxs = inputs_to_parameters!(state, io)
     sys = alias_elimination!(state)
     # TODO: avoid construct `TearingState` again.
     state = TearingState(sys)
@@ -984,7 +984,7 @@ end
 
 function io_preprocessing(sys::AbstractSystem, inputs,
                           outputs; simplify = false, kwargs...)
-    sys, input_idxs = structural_simplify(sys, (inputs, outputs); simplify, kwargs...)
+    sys, input_idxs = structural_simplify(sys, (; inputs, outputs); simplify, kwargs...)
 
     eqs = equations(sys)
     alg_start_idx = findfirst(!isdiffeq, eqs)
@@ -1201,6 +1201,8 @@ function linearize(sys, lin_fun; t = 0.0, op = Dict(), allow_input_derivatives =
     nz = size(f_z, 2)
     ny = size(h_x, 1)
 
+    D = h_u
+
     if isempty(g_z)
         A = f_x
         B = f_u
@@ -1216,20 +1218,20 @@ function linearize(sys, lin_fun; t = 0.0, op = Dict(), allow_input_derivatives =
         A = [f_x f_z
              gzgx*f_x gzgx*f_z]
         B = [f_u
-             zeros(nz, nu)]
+             gzgx * f_u] # The cited paper has zeros in the bottom block, see derivation in https://github.com/SciML/ModelingToolkit.jl/pull/1691 for the correct formula
+
         C = [h_x h_z]
         Bs = -(gz \ g_u) # This equation differ from the cited paper, the paper is likely wrong since their equaiton leads to a dimension mismatch.
         if !iszero(Bs)
             if !allow_input_derivatives
                 der_inds = findall(vec(any(!=(0), Bs, dims = 1)))
-                error("Input derivatives appeared in expressions (-g_z\\g_u != 0), the following inputs appeared differentiated: $(inputs(sys)[der_inds]). Call `linear_staespace` with keyword argument `allow_input_derivatives = true` to allow this and have the returned `B` matrix be of double width ($(2nu)), where the last $nu inputs are the derivatives of the first $nu inputs.")
+                error("Input derivatives appeared in expressions (-g_z\\g_u != 0), the following inputs appeared differentiated: $(inputs(sys)[der_inds]). Call `linear_statespace` with keyword argument `allow_input_derivatives = true` to allow this and have the returned `B` matrix be of double width ($(2nu)), where the last $nu inputs are the derivatives of the first $nu inputs.")
             end
-            B = [B Bs]
+            B = [B [zeros(nx, nu); Bs]]
+            D = [D zeros(ny, nu)]
         end
     end
 
-    D = h_u
-
     (; A, B, C, D)
 end
 

src/systems/diffeqs/odesystem.jl

@@ -316,7 +316,7 @@ function build_explicit_observed_function(sys, ts;
                     subs[s] = s′
                     continue
                 end
-                throw(ArgumentError("$s is either an observed nor a state variable."))
+                throw(ArgumentError("$s is neither an observed nor a state variable."))
             end
             continue
         end

src/systems/systemstructure.jl

@@ -403,11 +403,93 @@ function linear_subsys_adjmat(state::TransformationState)
                             linear_equations, eadj, cadj)
 end
 
+using .BipartiteGraphs: Label, BipartiteAdjacencyList
+struct SystemStructurePrintMatrix <:
+       AbstractMatrix{Union{Label, Int, BipartiteAdjacencyList}}
+    bpg::BipartiteGraph
+    highlight_graph::BipartiteGraph
+    var_to_diff::DiffGraph
+    eq_to_diff::DiffGraph
+    var_eq_matching::Union{Matching, Nothing}
+end
+Base.size(bgpm::SystemStructurePrintMatrix) = (max(nsrcs(bgpm.bpg), ndsts(bgpm.bpg)) + 1, 5)
+function compute_diff_label(diff_graph, i)
+    di = i - 1 <= length(diff_graph) ? diff_graph[i - 1] : nothing
+    ii = i - 1 <= length(invview(diff_graph)) ? invview(diff_graph)[i - 1] : nothing
+    return Label(string(di === nothing ? "" : string(di, '↓'),
+                        di !== nothing && ii !== nothing ? " " : "",
+                        ii === nothing ? "" : string(ii, '↑')))
+end
+function Base.getindex(bgpm::SystemStructurePrintMatrix, i::Integer, j::Integer)
+    checkbounds(bgpm, i, j)
+    if i <= 1
+        return (Label.(("#", "∂ₜ", "eq", "∂ₜ", "v")))[j]
+    elseif j == 2
+        return compute_diff_label(bgpm.eq_to_diff, i)
+    elseif j == 4
+        return compute_diff_label(bgpm.var_to_diff, i)
+    elseif j == 1
+        return i - 1
+    elseif j == 3
+        return BipartiteAdjacencyList(i - 1 <= nsrcs(bgpm.bpg) ?
+                                      𝑠neighbors(bgpm.bpg, i - 1) : nothing,
+                                      bgpm.highlight_graph !== nothing &&
+                                      i - 1 <= nsrcs(bgpm.highlight_graph) ?
+                                      Set(𝑠neighbors(bgpm.highlight_graph, i - 1)) :
+                                      nothing,
+                                      bgpm.var_eq_matching !== nothing &&
+                                      (i - 1 <= length(invview(bgpm.var_eq_matching))) ?
+                                      invview(bgpm.var_eq_matching)[i - 1] : unassigned)
+    elseif j == 5
+        return BipartiteAdjacencyList(i - 1 <= ndsts(bgpm.bpg) ?
+                                      𝑑neighbors(bgpm.bpg, i - 1) : nothing,
+                                      bgpm.highlight_graph !== nothing &&
+                                      i - 1 <= ndsts(bgpm.highlight_graph) ?
+                                      Set(𝑑neighbors(bgpm.highlight_graph, i - 1)) :
+                                      nothing,
+                                      bgpm.var_eq_matching !== nothing &&
+                                      (i - 1 <= length(bgpm.var_eq_matching)) ?
+                                      bgpm.var_eq_matching[i - 1] : unassigned)
+    else
+        @assert false
+    end
+end
+
 function Base.show(io::IO, mime::MIME"text/plain", s::SystemStructure)
-    @unpack graph = s
-    S = incidence_matrix(graph, Num(Sym{Real}(:×)))
-    print(io, "Incidence matrix:")
-    show(io, mime, S)
+    @unpack graph, solvable_graph, var_to_diff, eq_to_diff = s
+    if !get(io, :limit, true) || !get(io, :mtk_limit, true)
+        print(io, "SystemStructure with ", length(graph.fadjlist), " equations and ",
+              isa(graph.badjlist, Int) ? graph.badjlist : length(graph.badjlist),
+              " variables\n")
+        Base.print_matrix(io,
+                          SystemStructurePrintMatrix(complete(graph),
+                                                     complete(solvable_graph),
+                                                     complete(var_to_diff),
+                                                     complete(eq_to_diff), nothing))
+    else
+        S = incidence_matrix(graph, Num(Sym{Real}(:×)))
+        print(io, "Incidence matrix:")
+        show(io, mime, S)
+    end
+end
+
+struct MatchedSystemStructure
+    structure::SystemStructure
+    var_eq_matching::Matching
+end
+
+function Base.show(io::IO, mime::MIME"text/plain", ms::MatchedSystemStructure)
+    s = ms.structure
+    @unpack graph, solvable_graph, var_to_diff, eq_to_diff = s
+    print(io, "Matched SystemStructure with ", length(graph.fadjlist), " equations and ",
+          isa(graph.badjlist, Int) ? graph.badjlist : length(graph.badjlist),
+          " variables\n")
+    Base.print_matrix(io,
+                      SystemStructurePrintMatrix(complete(graph),
+                                                 complete(solvable_graph),
+                                                 complete(var_to_diff),
+                                                 complete(eq_to_diff),
+                                                 complete(ms.var_eq_matching, nsrcs(graph))))
 end
 
 end # module

test/structural_transformation/tearing.jl

@@ -30,6 +30,12 @@ if VERSION >= v"1.6"
     @test occursin("Incidence matrix:", prt)
     @test occursin("×", prt)
     @test occursin("⋅", prt)
+
+    buff = IOBuffer()
+    io = IOContext(buff, :mtk_limit => false)
+    show(io, MIME"text/plain"(), state.structure)
+    prt = String(take!(buff))
+    @test occursin("SystemStructure", prt)
 end
 
 # u1 = f1(u5)