refactor: move full_equations to ModelingToolkit, don't error on singular systems

Author: AayushSabharwal

src/structural_transformation/symbolics_tearing.jl

@@ -99,48 +99,6 @@ function eq_derivative!(ts::TearingState, ieq::Int; kwargs...)
     return eq_diff
 end
 
-function tearing_sub(expr, dict, s)
-    expr = Symbolics.fixpoint_sub(expr, dict; operator = ModelingToolkit.Initial)
-    s ? simplify(expr) : expr
-end
-
-function tearing_substitute_expr(sys::AbstractSystem, expr; simplify = false)
-    empty_substitutions(sys) && return expr
-    substitutions = get_substitutions(sys)
-    return tearing_sub(expr, substitutions, simplify)
-end
-
-"""
-$(TYPEDSIGNATURES)
-
-Like `equations(sys)`, but includes substitutions done by the tearing process.
-These equations matches generated numerical code.
-
-See also [`equations`](@ref) and [`ModelingToolkit.get_eqs`](@ref).
-"""
-function full_equations(sys::AbstractSystem; simplify = false)
-    isempty(observed(sys)) && return equations(sys)
-    subs = Dict([eq.lhs => eq.rhs for eq in observed(sys)])
-    neweqs = map(equations(sys)) do eq
-        if iscall(eq.lhs) && operation(eq.lhs) isa Union{Shift, Differential}
-            return tearing_sub(eq.lhs, subs, simplify) ~ tearing_sub(eq.rhs, subs,
-                simplify)
-        else
-            if !(eq.lhs isa Number && eq.lhs == 0)
-                eq = 0 ~ eq.rhs - eq.lhs
-            end
-            rhs = tearing_sub(eq.rhs, subs, simplify)
-            if rhs isa Symbolic
-                return 0 ~ rhs
-            else # a number
-                error("tearing failed because the system is singular")
-            end
-        end
-        eq
-    end
-    return neweqs
-end
-
 function tearing_substitution(sys::AbstractSystem; kwargs...)
     neweqs = full_equations(sys::AbstractSystem; kwargs...)
     @set! sys.eqs = neweqs

src/systems/abstractsystem.jl

@@ -1679,6 +1679,57 @@ function equations_toplevel(sys::AbstractSystem)
     return get_eqs(sys)
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Recursively substitute `dict` into `expr`. Use `Symbolics.simplify` on the expression
+if `simplify == true`.
+"""
+function substitute_and_simplify(expr, dict::AbstractDict, simplify::Bool)
+    expr = Symbolics.fixpoint_sub(expr, dict; operator = ModelingToolkit.Initial)
+    simplify ? Symbolics.simplify(expr) : expr
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Recursively substitute the observed equations of `sys` into `expr`. If `simplify`, call
+`Symbolics.simplify` on the resultant expression.
+"""
+function substitute_observed(sys::AbstractSystem, expr; simplify = false)
+    empty_substitutions(sys) && return expr
+    substitutions = get_substitutions(sys)
+    return substitute_and_simplify(expr, substitutions, simplify)
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Like `equations(sys)`, but also substitutes the observed equations eliminated from the
+equations during `mtkcompile`. These equations matches generated numerical code.
+
+See also [`equations`](@ref) and [`ModelingToolkit.get_eqs`](@ref).
+"""
+function full_equations(sys::AbstractSystem; simplify = false)
+    empty_substitutions(sys) && return equations(sys)
+    subs = get_substitutions(sys)
+    neweqs = map(equations(sys)) do eq
+        if iscall(eq.lhs) && operation(eq.lhs) isa Union{Shift, Differential}
+            return substitute_and_simplify(eq.lhs, subs, simplify) ~ substitute_and_simplify(
+                eq.rhs, subs,
+                simplify)
+        else
+            if !_iszero(eq.lhs)
+                eq = 0 ~ eq.rhs - eq.lhs
+            end
+            rhs = substitute_and_simplify(eq.rhs, subs, simplify)
+            return 0 ~ rhs
+        end
+        eq
+    end
+    return neweqs
+end
+
 """
     $(TYPEDSIGNATURES)