Update alias_elimination tests

Author: YingboMa

src/systems/reduction.jl

@@ -33,58 +33,6 @@ isvar(s::Sym; param=false) = param ? true : !isparameter(s)
 isvar(s::Term; param=false) = isvar(s.op; param=param)
 isvar(s::Any;param=false) = false
 
-function filterexpr(f, s)
-    vs = []
-    Rewriters.Prewalk(Rewriters.Chain([@rule((~x::f) => push!(vs, ~x))]))(s)
-    vs
-end
-
-function make_lhs_0(eq)
-    if eq.lhs isa Number && iszero(eq.lhs)
-        return eq
-    else
-        0 ~ eq.lhs - eq.rhs
-    end
-end
-
-function alias_elimination(sys::ODESystem)
-    eqs = vcat(equations(sys), observed(sys))
-
-    # make all algebraic equations have 0 on LHS
-    eqs = map(eqs) do eq
-        if eq.lhs isa Term && eq.lhs.op isa Differential
-            eq
-        else
-            make_lhs_0(eq)
-        end
-    end
-
-    newstates = map(eqs) do eq
-        if eq.lhs isa Term && eq.lhs.op isa Differential
-            filterexpr(isvar, eq.lhs)
-        else
-            []
-        end
-    end |> Iterators.flatten |> collect |> unique
-
-
-    all_vars = map(eqs) do eq
-        filterexpr(isvar, eq.rhs)
-    end |> Iterators.flatten |> collect |> unique
-
-    alg_idxs = findall(x->!(x.lhs isa Term) && iszero(x.lhs), eqs)
-
-    eliminate = setdiff(all_vars, newstates)
-
-    outputs = solve_for(eqs[alg_idxs], eliminate)
-
-    diffeqs = eqs[setdiff(1:length(eqs), alg_idxs)]
-
-    diffeqs′ = substitute_aliases(diffeqs, Dict(eliminate .=> outputs))
-
-    ODESystem(diffeqs′, sys.iv, newstates, parameters(sys), observed=eliminate .~ outputs)
-end
-
 function get_α_x(αx)
     if isvar(αx, param=true)
         return 1, αx
@@ -101,7 +49,7 @@ function get_α_x(αx)
     end
 end
 
-function alias_elimination2(sys)
+function alias_elimination(sys)
     eqs = vcat(equations(sys), observed(sys))
 
     subs = Pair[]

test/reduction.jl

@@ -1,5 +1,4 @@
 using ModelingToolkit, OrdinaryDiffEq, Test
-using ModelingToolkit: alias_elimination2
 
 @parameters t σ ρ β
 @variables x(t) y(t) z(t) a(t) u(t) F(t)
@@ -15,7 +14,7 @@ eqs = [D(x) ~ σ*(y-x),
 
 lorenz1 = ODESystem(eqs,t,[x,y,z,a],[σ,ρ,β],name=:lorenz1)
 
-lorenz1_aliased = alias_elimination2(lorenz1)
+lorenz1_aliased = alias_elimination(lorenz1)
 @test length(equations(lorenz1_aliased)) == 3
 @test length(states(lorenz1_aliased)) == 3
 
@@ -47,7 +46,7 @@ lorenz2 = ODESystem(eqs2,pins=[F],observed=aliases2,name=:lorenz2)
 connections = [lorenz1.F ~ lorenz2.u,
                lorenz2.F ~ lorenz1.u]
 
-connected = ODESystem([lorenz2.y ~ a + lorenz1.x ],t,[a],[],observed=connections,systems=[lorenz1,lorenz2])
+connected = ODESystem([lorenz2.y ~ a + lorenz1.x],t,[a],[],observed=connections,systems=[lorenz1,lorenz2])
 
 # Reduced Unflattened System
 #=
@@ -62,7 +61,7 @@ connected = ODESystem(Equation[],t,[],[],observed=connections2,systems=[lorenz1,
 
 flattened_system = ModelingToolkit.flatten(connected)
 
-aliased_flattened_system = alias_elimination2(flattened_system)
+aliased_flattened_system = alias_elimination(flattened_system)
 
 @test isequal(states(aliased_flattened_system), [
         lorenz1.x
@@ -83,21 +82,24 @@ aliased_flattened_system = alias_elimination2(flattened_system)
         lorenz2.β
        ]) |> isempty
 
-test_equal.(equations(aliased_flattened_system), [
-       D(lorenz1.x) ~ lorenz1.σ*(lorenz1.y-lorenz1.x) + lorenz2.x - lorenz2.y - lorenz2.z,
-       D(lorenz1.y) ~ lorenz1.x*(lorenz1.ρ-lorenz1.z)-(lorenz1.x + lorenz1.y - lorenz1.z),
-       D(lorenz1.z) ~ lorenz1.x*lorenz1.y - lorenz1.β*lorenz1.z,
-       D(lorenz2.x) ~ lorenz1.x + lorenz1.y - lorenz1.z,
-       D(lorenz2.y) ~ lorenz2.x*(lorenz2.ρ-lorenz2.z)-lorenz2.x,
-       D(lorenz2.z) ~ lorenz2.x*lorenz2.y - lorenz2.β*lorenz2.z])
-
-test_equal.(observed(aliased_flattened_system), [
-       lorenz1.F ~ lorenz2.x + -1 * (lorenz2.y + lorenz2.z),
-       lorenz1.u ~ lorenz1.x + lorenz1.y + -1 * lorenz1.z,
-       lorenz2.F ~ lorenz1.x + lorenz1.y + -1 * lorenz1.z,
-       a ~ lorenz2.y + -1 * lorenz1.x,
-       lorenz2.u ~ lorenz2.x + -1 * (lorenz2.y + lorenz2.z),
-])
+reduced_eqs = [
+               D(lorenz1.x) ~ lorenz1.σ*(lorenz1.y-lorenz1.x) + lorenz2.x - (a + lorenz1.x) - lorenz2.z,
+               D(lorenz1.y) ~ lorenz1.x*(lorenz1.ρ-lorenz1.z)-(lorenz1.x + lorenz1.y - lorenz1.z),
+               D(lorenz1.z) ~ lorenz1.x*lorenz1.y - lorenz1.β*lorenz1.z,
+               D(lorenz2.x) ~ lorenz1.x + lorenz1.y - lorenz1.z,
+               D(lorenz2.y) ~ lorenz2.x*(lorenz2.ρ-lorenz2.z)-lorenz2.x,
+               D(lorenz2.z) ~ lorenz2.x*(a + lorenz1.x) - lorenz2.β*lorenz2.z
+              ]
+test_equal.(equations(aliased_flattened_system), reduced_eqs)
+
+observed_eqs = [
+                lorenz2.y ~ a + lorenz1.x,
+                lorenz1.F ~ lorenz2.u,
+                lorenz2.F ~ lorenz1.u,
+                lorenz1.u ~ lorenz1.x + lorenz1.y - lorenz1.z,
+                lorenz2.u ~ lorenz2.x - lorenz2.y - lorenz2.z,
+               ]
+test_equal.(observed(aliased_flattened_system), observed_eqs)
 
 
 # issue #578
@@ -110,7 +112,7 @@ let
         x ~ y
     ];
     sys = ODESystem(eqs, t, [x], []);
-    asys = alias_elimination2(ModelingToolkit.flatten(sys))
+    asys = alias_elimination(ModelingToolkit.flatten(sys))
 
     test_equal.(asys.eqs, [D(x) ~ 2x])
     test_equal.(asys.observed, [y ~ x])