fixpoint the substitution of variables

Author: shashi

src/systems/reduction.jl

@@ -3,14 +3,24 @@ export alias_elimination
 function flatten(sys::ODESystem)
     ODESystem(equations(sys),
               independent_variable(sys),
-              states(sys),
-              parameters(sys),
               observed=observed(sys))
 end
 
 
+using SymbolicUtils: Rewriters
+
+function fixpoint_sub(x, dict)
+    y = substitute(x, dict)
+    while !isequal(x, y)
+        y = x
+        x = substitute(y, dict)
+    end
+
+    return x
+end
+
 function substitute_aliases(diffeqs, outputs)
-    lhss(diffeqs) .~ substitute.(rhss(diffeqs), (Dict(lhss(outputs) .=> rhss(outputs)),))
+    lhss(diffeqs) .~ fixpoint_sub.(rhss(diffeqs), (Dict(lhss(outputs) .=> rhss(outputs)),))
 end
 
 function make_lhs_0(eq)
@@ -32,11 +42,17 @@ function alias_elimination(sys::ODESystem)
             []
         end
     end |> Iterators.flatten |> collect |> unique
+
+    all_vars = map(eqs) do eq
+        filter(x->!isparameter(x.op), get_variables(eq.rhs))
+    end |> Iterators.flatten |> collect |> unique
+
     newstates = convert.(Variable, new_stateops)
 
     alg_idxs = findall(x->x.lhs isa Constant && iszero(x.lhs), eqs)
 
-    eliminate = setdiff(states(sys), newstates)
+    eliminate = setdiff(convert.(Variable, all_vars), newstates)
+
     outputs = solve_for(eqs[alg_idxs], map(x->x(sys.iv()), eliminate))
 
     diffeqs = eqs[setdiff(1:length(eqs), alg_idxs)]

test/reduction.jl

@@ -1,9 +1,11 @@
 using ModelingToolkit, OrdinaryDiffEq, Test
 
-@parameters t σ ρ β F(t)
-@variables x(t) y(t) z(t) a(t) u(t)
+@parameters t σ ρ β
+@variables x(t) y(t) z(t) a(t) u(t) F(t)
 @derivatives D'~t
 
+test_equal(a, b) = @test isequal(simplify(a, polynorm=true), simplify(b, polynorm=true))
+
 eqs = [D(x) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,
        D(z) ~ a*y - β*z,
@@ -78,16 +80,16 @@ aliased_flattened_system = alias_elimination(flattened_system)
         lorenz2.β
        ])) |> isempty
 
-#=
-equations(reduced_flattened_system) == [
+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*(ρ-z)-lorenz1.x - lorenz1.y + lorenz1.z,
-       D(lorenz1.z) ~ lorenz1.x*lorenz1.y - lorenz1.β*lorenz1.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]
+       D(lorenz2.z) ~ lorenz2.x*lorenz2.y - lorenz2.β*lorenz2.z])
 
-observed(reduced_flattened_system) == [
+#=
+observed(aliased_flattened_system) == [
         lorenz1.F ~ lorenz2.u
         lorenz2.F ~ lorenz1.u
         a ~ -lorenz1.x + lorenz2.y