Merge pull request #1208 from SciML/myb/opt

Don't reduce inputs

Author: YingboMa

src/ModelingToolkit.jl

@@ -52,7 +52,7 @@ import Symbolics: rename, get_variables!, _solve, hessian_sparsity,
                   ParallelForm, SerialForm, MultithreadedForm, build_function,
                   unflatten_long_ops, rhss, lhss, prettify_expr, gradient,
                   jacobian, hessian, derivative, sparsejacobian, sparsehessian,
-                  substituter, scalarize
+                  substituter, scalarize, getparent
 
 import DiffEqBase: @add_kwonly
 

src/structural_transformation/StructuralTransformations.jl

@@ -16,7 +16,8 @@ using ModelingToolkit
 using ModelingToolkit: ODESystem, AbstractSystem,var_from_nested_derivative, Differential,
                        states, equations, vars, Symbolic, diff2term, value,
                        operation, arguments, Sym, Term, simplify, solve_for,
-                       isdiffeq, isdifferential, get_structure, get_iv, independent_variables,
+                       isdiffeq, isdifferential, isinput,
+                       get_structure, get_iv, independent_variables,
                        get_structure, defaults, InvalidSystemException,
                        ExtraEquationsSystemException,
                        ExtraVariablesSystemException,

src/structural_transformation/utils.jl

@@ -276,7 +276,9 @@ function find_solvables!(sys)
         term = value(eq.rhs - eq.lhs)
         for j in 𝑠neighbors(graph, i)
             isalgvar(s, j) || continue
-            a, b, islinear = linear_expansion(term, fullvars[j])
+            var = fullvars[j]
+            isinput(var) && continue
+            a, b, islinear = linear_expansion(term, var)
             a = unwrap(a)
             if islinear && (!(a isa Symbolic) && a isa Number && a != 0)
                 add_edge!(solvable_graph, i, j)

src/systems/systemstructure.jl

@@ -5,7 +5,8 @@ using Symbolics: linear_expansion, unwrap
 using SymbolicUtils: istree, operation, arguments, Symbolic
 using ..ModelingToolkit
 import ..ModelingToolkit: isdiffeq, var_from_nested_derivative, vars!, flatten,
-    value, InvalidSystemException, isdifferential, _iszero, isparameter, independent_variables
+    value, InvalidSystemException, isdifferential, _iszero, isparameter,
+    independent_variables, isinput
 using ..BipartiteGraphs
 using LightGraphs
 using UnPack
@@ -230,7 +231,7 @@ function find_linear_equations(sys)
             var = fullvars[j]
             a, b, islinear = linear_expansion(term, var)
             a = unwrap(a)
-            if islinear && !(a isa Symbolic) && a isa Number
+            if islinear && !(a isa Symbolic) && a isa Number && !isinput(var)
                 if a == 1 || a == -1
                     a = convert(Integer, a)
                     linear_term += a * var

src/utils.jl

@@ -225,7 +225,12 @@ isdiffeq(eq) = isdifferential(eq.lhs)
 isdifference(expr) = istree(expr) && operation(expr) isa Difference
 isdifferenceeq(eq) = isdifference(eq.lhs)
 
-isvariable(x) = x isa Symbolic && hasmetadata(x, VariableSource)
+function isvariable(x)
+    x isa Symbolic || return false
+    p = getparent(x, nothing)
+    p === nothing || (x = p)
+    hasmetadata(x, VariableSource)
+end
 
 vars(x::Sym; op=Differential) = Set([x])
 vars(exprs::Symbolic; op=Differential) = vars([exprs]; op=op)

src/variables.jl

@@ -11,6 +11,15 @@ Symbolics.option_to_metadata_type(::Val{:description}) = VariableDescriptionType
 Symbolics.option_to_metadata_type(::Val{:input}) = VariableInput
 Symbolics.option_to_metadata_type(::Val{:output}) = VariableOutput
 
+function isvarkind(m, x)
+    p = getparent(x, nothing)
+    p === nothing || (x = p)
+    getmetadata(x, m, false)
+end
+
+isinput(x) = isvarkind(VariableInput, x)
+isoutput(x) = isvarkind(VariableOutput, x)
+
 """
 $(SIGNATURES)
 

test/odesystem.jl

@@ -399,6 +399,7 @@ eqs = [
 @test isequal(@nonamespace(sys.y), unwrap(y))
 @test isequal(@nonamespace(sys.p), unwrap(p))
 @test_nowarn sys.x, sys.y, sys.p
+@test ModelingToolkit.isvariable(Symbolics.unwrap(x[1]))
 
 # Mixed Difference Differential equations
 @parameters t a b c d

test/reduction.jl

@@ -281,3 +281,20 @@ dt = ModelingToolkit.value(ModelingToolkit.derivative(eq.rhs, sin(10t)))
 @test dv25 ≈ 0.3
 @test ddv25 == 0.005
 @test dt == -0.1
+
+# Don't reduce inputs
+@parameters t σ ρ β
+@variables x(t) y(t) z(t) [input=true] a(t) u(t) F(t)
+D = Differential(t)
+
+eqs = [
+       D(x) ~ σ*(y-x)
+       D(y) ~ x*(ρ-z)-y + β
+       0 ~ z - x + y
+       0 ~ a + z
+       u ~ z + a
+      ]
+
+lorenz1 = ODESystem(eqs,t,name=:lorenz1)
+lorenz1_reduced = structural_simplify(lorenz1)
+@test z in Set(states(lorenz1_reduced))