isisomorphic and convert_system fix

Author: anandijain

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"

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,53 @@ 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
+(f::Symbolic{<:FnType})(args...; kwargs...) = Term{promote_symtype(f, symtype.(args)...)}(f, [args...]; kwargs...)
+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
+            elseif !isequal(eq, eq2) && j == length(eqs2)
+            end
+        end
+    end
+    eqpairs
+end
+
+function isisomorphic(sys1::AbstractODESystem, sys2::AbstractODESystem; verbose=false)
+    sys1 = flatten(sys1)
+    sys2 = flatten(sys2)
+
+    iv1, iv2 = independent_variable(sys1), independent_variable(sys2) # not needed
+    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 = operation(s)(t; metadata=SymbolicUtils.metadata(s))
             newsts[i] = ns
             varmap[s] = ns
         else
@@ -377,7 +377,9 @@ function convert_system(::Type{<:ODESystem}, sys, t; name=nameof(sys))
         end
     end
     sub = Base.Fix2(substitute, varmap)
+    iv = independent_variable(sys)
+    sub.x[iv] = t # otherwise the Differentials aren't fixed
     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,45 @@ 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
+    ]
+
+    eqs3 = copy(eqs2)
+    
+    # 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)
+    
+    @named sys1 = ODESystem(eqs)
+    @named sys2 = ODESystem(eqs2)
+    @named sys3 = ODESystem(eqs3, t2)
+    ssys3 = structural_simplify(sys3)
+
+    @test ModelingToolkit.isisomorphic(sys1, sys2)
+    @test !ModelingToolkit.isisomorphic(sys1, sys3)
+    @test ModelingToolkit.isisomorphic(sys1, ssys3) # I don't call structural_simplify in isisomorphic
+
+    # 1281
+    iv2 = independent_variable(sys2)
+    @test isequal(independent_variable(convert_system(ODESystem, sys1, iv2)), iv2)
+end