Merge pull request #1538 from SciML/myb/option

Add allow_symbolic and allow_parameter as options in structural_simplify

Author: YingboMa

src/structural_transformation/pantelides.jl

@@ -139,7 +139,6 @@ instead, which calls this function internally.
 """
 function dae_index_lowering(sys::ODESystem; kwargs...)
     state = TearingState(sys)
-    find_solvables!(state)
     var_eq_matching = pantelides!(state; kwargs...)
     return invalidate_cache!(pantelides_reassemble(state, var_eq_matching))
 end

src/structural_transformation/symbolics_tearing.jl

@@ -33,7 +33,7 @@ function var_derivative!(ts::TearingState{ODESystem}, v::Int)
     sys = ts.sys
     s = ts.structure
     D = Differential(get_iv(sys))
-    add_vertex!(s.solvable_graph, DST)
+    s.solvable_graph === nothing || add_vertex!(s.solvable_graph, DST)
     push!(ts.fullvars, D(ts.fullvars[v]))
 end
 
@@ -43,7 +43,7 @@ function eq_derivative!(ts::TearingState{ODESystem}, ieq::Int)
     D = Differential(get_iv(sys))
     eq = equations(ts)[ieq]
     eq = ModelingToolkit.expand_derivatives(0 ~ D(eq.rhs - eq.lhs))
-    add_vertex!(s.solvable_graph, SRC)
+    s.solvable_graph === nothing || add_vertex!(s.solvable_graph, SRC)
     push!(equations(ts), eq)
     # Analyze the new equation and update the graph/solvable_graph
     # First, copy the previous incidence and add the derivative terms.
@@ -54,7 +54,7 @@ function eq_derivative!(ts::TearingState{ODESystem}, ieq::Int)
         add_edge!(s.graph, eq_diff, var)
         add_edge!(s.graph, eq_diff, s.var_to_diff[var])
     end
-    find_eq_solvables!(ts, eq_diff; may_be_zero=true, allow_symbolic=true)
+    s.solvable_graph === nothing || find_eq_solvables!(ts, eq_diff; may_be_zero=true, allow_symbolic=true)
 end
 
 function tearing_sub(expr, dict, s)
@@ -242,8 +242,8 @@ function tearing_reassemble(state::TearingState, var_eq_matching; simplify=false
     return invalidate_cache!(sys)
 end
 
-function tearing(state::TearingState)
-    state.structure.solvable_graph === nothing && find_solvables!(state)
+function tearing(state::TearingState; kwargs...)
+    state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
     complete!(state.structure)
     @unpack graph, solvable_graph = state.structure
     algvars = BitSet(findall(v->isalgvar(state.structure, v), 1:ndsts(graph)))

src/structural_transformation/utils.jl

@@ -197,13 +197,13 @@ function find_eq_solvables!(state::TearingState, ieq; may_be_zero=false, allow_s
     end
 end
 
-function find_solvables!(state::TearingState; allow_symbolic=false)
+function find_solvables!(state::TearingState; kwargs...)
     @assert state.structure.solvable_graph === nothing
     eqs = equations(state)
     graph = state.structure.graph
     state.structure.solvable_graph = BipartiteGraph(nsrcs(graph), ndsts(graph))
     for ieq in 1:length(eqs)
-        find_eq_solvables!(state, ieq; allow_symbolic)
+        find_eq_solvables!(state, ieq; kwargs...)
     end
     return nothing
 end

src/systems/abstractsystem.jl

@@ -919,9 +919,11 @@ $(SIGNATURES)
 
 Structurally simplify algebraic equations in a system and compute the
 topological sort of the observed equations. When `simplify=true`, the `simplify`
-function will be applied during the tearing process.
+function will be applied during the tearing process. It also takes kwargs
+`allow_symbolic=false` and `allow_parameter=true` which limits the coefficient
+types during tearing.
 """
-function structural_simplify(sys::AbstractSystem; simplify=false)
+function structural_simplify(sys::AbstractSystem; simplify=false, kwargs...)
     sys = expand_connections(sys)
     sys = alias_elimination(sys)
     state = TearingState(sys)
@@ -930,7 +932,7 @@ function structural_simplify(sys::AbstractSystem; simplify=false)
         sys = dae_index_lowering(ode_order_lowering(sys))
     end
     state = TearingState(sys)
-    find_solvables!(state)
+    find_solvables!(state; kwargs...)
     sys = tearing_reassemble(state, tearing(state), simplify=simplify)
     fullstates = [map(eq->eq.lhs, observed(sys)); states(sys)]
     @set! sys.observed = topsort_equations(observed(sys), fullstates)

test/components.jl

@@ -4,10 +4,10 @@ using ModelingToolkit.BipartiteGraphs
 using ModelingToolkit.StructuralTransformations
 
 function check_contract(sys)
+    graph = ModelingToolkit.get_tearing_state(sys).structure.graph
     sys = tearing_substitution(sys)
     state = TearingState(sys)
     fullvars = state.fullvars
-    graph = state.structure.graph
 
     eqs = equations(sys)
     var2idx = Dict(enumerate(fullvars))
@@ -30,6 +30,7 @@ end
 
 include("../examples/rc_model.jl")
 
+@test length(equations(structural_simplify(rc_model, allow_parameter=false))) > 1
 sys = structural_simplify(rc_model)
 check_contract(sys)
 @test !isempty(ModelingToolkit.defaults(sys))