Merge pull request #1543 from anandijain/aj/iso2

isisomorphic and convert_system fix

Author: YingboMa

Project.toml

@@ -6,6 +6,7 @@ version = "8.8.1"
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
@@ -45,6 +46,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 AbstractTrees = "0.3"
 ArrayInterface = "3.1.39, 4, 5"
+Combinatorics = "1"
 ConstructionBase = "1"
 DataStructures = "0.17, 0.18"
 DiffEqBase = "6.81.0"

src/ModelingToolkit.jl

@@ -23,6 +23,7 @@ import MacroTools: splitdef, combinedef, postwalk, striplines
 import Libdl
 using DocStringExtensions
 using Base: RefValue
+using Combinatorics
 import IfElse
 
 import Distributions

src/systems/diffeqs/abstractodesystem.jl

@@ -900,3 +900,51 @@ end
 function SteadyStateProblemExpr(sys::AbstractODESystem, args...; kwargs...)
     SteadyStateProblemExpr{true}(sys, args...; kwargs...)
 end
+
+# SU 442 and Symbolics 588 copied over for CI, will be removed
+function SymbolicUtils.substitute(op::Symbolics.Differential, dict; kwargs...)
+    @set! op.x = substitute(op.x, dict; kwargs...)
+end
+
+function _match_eqs(eqs1, eqs2)
+    eqpairs = Pair[]
+    for (i, eq) in enumerate(eqs1)
+        for (j, eq2) in enumerate(eqs2)
+            if isequal(eq, eq2)
+                push!(eqpairs, i => j)
+                break
+            end
+        end
+    end
+    eqpairs
+end
+
+function isisomorphic(sys1::AbstractODESystem, sys2::AbstractODESystem)
+    sys1 = flatten(sys1)
+    sys2 = flatten(sys2)
+
+    iv2 = only(independent_variables(sys2))
+    sys1 = convert_system(ODESystem, sys1, iv2)
+    s1, s2 = states(sys1), states(sys2)
+    p1, p2 = parameters(sys1), parameters(sys2)
+
+    (length(s1) != length(s2)) || (length(p1) != length(p2)) && return false
+
+    eqs1 = equations(sys1)
+    eqs2 = equations(sys2)
+
+    pps = permutations(p2)
+    psts = permutations(s2)
+    orig = [p1; s1]
+    perms = ([x; y] for x in pps for y in psts)
+
+    for perm in perms
+        rules = Dict(orig .=> perm)
+        neweqs1 = substitute(eqs1, rules)
+        eqpairs = _match_eqs(neweqs1, eqs2)
+        if length(eqpairs) == length(eqs1)
+            return true
+        end
+    end
+    return false
+end

src/systems/diffeqs/odesystem.jl

@@ -367,7 +367,7 @@ function convert_system(::Type{<:ODESystem}, sys, t; name=nameof(sys))
                 newsts[i] = s
                 continue
             end
-            ns = operation(s)(t)
+            ns = similarterm(s, operation(s), (t,); metadata=SymbolicUtils.metadata(s))
             newsts[i] = ns
             varmap[s] = ns
         else
@@ -377,7 +377,11 @@ function convert_system(::Type{<:ODESystem}, sys, t; name=nameof(sys))
         end
     end
     sub = Base.Fix2(substitute, varmap)
+    if sys isa AbstractODESystem
+        iv = only(independent_variables(sys))
+        sub.x[iv] = t # otherwise the Differentials aren't fixed
+    end
     neweqs = map(sub, equations(sys))
     defs = Dict(sub(k) => sub(v) for (k, v) in defaults(sys))
-    return ODESystem(neweqs, t, newsts, parameters(sys); defaults=defs, name=name,checks=false)
+    return ODESystem(neweqs, t, newsts, parameters(sys); defaults=defs, name=name, checks=false)
 end

test/odesystem.jl

@@ -672,3 +672,51 @@ let
     @named sys = ODESystem(eqs, t)
     @test length(equations(structural_simplify(sys))) == 2
 end
+
+let 
+    eq_to_lhs(eq) = eq.lhs - eq.rhs ~ 0
+    eqs_to_lhs(eqs) = eq_to_lhs.(eqs)
+
+    @parameters σ = 10 ρ = 28 β = 8 / 3 sigma rho beta
+    @variables t t2 x(t) = 1 y(t) = 0 z(t) = 0 x2(t2) = 1 y2(t2) = 0 z2(t2) = 0 u[1:3](t2)
+
+    D = Differential(t)
+    D2 = Differential(t2)
+
+    eqs = [D(x) ~ σ * (y - x),
+        D(y) ~ x * (ρ - z) - y,
+        D(z) ~ x * y - β * z]
+
+    eqs2 = [
+        D2(y2) ~ x2 * (rho - z2) - y2,
+        D2(x2) ~ sigma * (y2 - x2),
+        D2(z2) ~ x2 * y2 - beta * z2
+    ]
+
+    # array u
+    eqs3 = [D2(u[1]) ~ sigma * (u[2] - u[1]),
+        D2(u[2]) ~ u[1] * (rho - u[3]) - u[2],
+        D2(u[3]) ~ u[1] * u[2] - beta * u[3]]
+    eqs3 = eqs_to_lhs(eqs3)
+    
+    eqs4 = [
+        D2(y2) ~ x2 * (rho - z2) - y2,
+        D2(x2) ~ sigma * (y2 - x2),
+        D2(z2) ~ y2 - beta * z2 # missing x2 term
+    ]
+
+    @named sys1 = ODESystem(eqs)
+    @named sys2 = ODESystem(eqs2)
+    @named sys3 = ODESystem(eqs3, t2)
+    ssys3 = structural_simplify(sys3)
+    @named sys4 = ODESystem(eqs4)
+
+    @test ModelingToolkit.isisomorphic(sys1, sys2)
+    @test !ModelingToolkit.isisomorphic(sys1, sys3)
+    @test ModelingToolkit.isisomorphic(sys1, ssys3) # I don't call structural_simplify in isisomorphic
+    @test !ModelingToolkit.isisomorphic(sys1, sys4) 
+
+    # 1281
+    iv2 = only(independent_variables(sys2))
+    @test isequal(only(independent_variables(convert_system(ODESystem, sys1, iv2))), iv2)
+end