Complete ModelingToolkit v9 and Julia v1.10 compatibility upgrade

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelOrderReduction"
 uuid = "207801d6-6cee-43a9-ad0c-f0c64933efa0"
-authors = ["Bowen S. Zhu <bszhu@fas.harvard.edu> and contributors"]
-version = "0.1.1"
+authors = ["Bowen S. Zhu <bowenzhu@mit.edu> and contributors"]
+version = "0.1.2"
 
 [deps]
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
@@ -10,18 +10,14 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 RandomizedLinAlg = "0448d7d9-159c-5637-8537-fd72090fea46"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
-Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
 
 [compat]
 DocStringExtensions = "0.8, 0.9"
 LinearAlgebra = "1"
-ModelingToolkit = "8.21"
+ModelingToolkit = "9"
 RandomizedLinAlg = "0.1"
 Setfield = "0.8, 1"
 SparseArrays = "1"
-SymbolicUtils = "0.19"
-Symbolics = "4.10"
 TSVD = "0.4"
-julia = "1.6"
+julia = "1.10"

docs/Project.toml

@@ -12,7 +12,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 DifferentialEquations = "7.6"
 Documenter = "1"
 LaTeXStrings = "1.3"
-MethodOfLines = "0.6, 0.7"
+MethodOfLines = "0.11"
 ModelOrderReduction = "0.1"
-ModelingToolkit = "8.33"
-Plots = "1.36"
+ModelingToolkit = "9"
+Plots = "1.40"

docs/pages.jl

@@ -1,5 +1,5 @@
 pages = [
     "Home" => "index.md",
     "Functions" => "functions.md",
-    "Tutorials" => Any["tutorials/deim_FitzHugh-Nagumo.md"]
+    "Tutorials" => Any["tutorials/deim_FitzHugh-Nagumo.md"],
 ]

src/ModelOrderReduction.jl

@@ -2,7 +2,6 @@ module ModelOrderReduction
 
 using DocStringExtensions
 
-using Symbolics
 using ModelingToolkit
 using LinearAlgebra
 

src/deim.jl

@@ -130,22 +130,22 @@ function deim(sys::ODESystem, snapshot::AbstractMatrix, pod_dim::Integer;
 
     iv = ModelingToolkit.get_iv(sys) # the single independent variable
     D = Differential(iv)
-    dvs = ModelingToolkit.get_states(sys) # dependent variables
+    dvs = ModelingToolkit.get_unknowns(sys) # dependent variables
 
     pod_reducer = POD(snapshot, pod_dim)
     reduce!(pod_reducer, TSVD())
     V = pod_reducer.rbasis # POD basis
 
     var_name = gensym(:ŷ)
     ŷ = (@variables $var_name(iv)[1:pod_dim])[1]
-    @set! sys.states = Symbolics.value.(Symbolics.scalarize(ŷ)) # new variables from POD
+    @set! sys.unknowns = Symbolics.value.(Symbolics.scalarize(ŷ)) # new variables from POD
     ModelingToolkit.get_var_to_name(sys)[Symbolics.getname(ŷ)] = Symbolics.unwrap(ŷ)
 
     deqs, eqs = get_deqs(sys) # split eqs into differential and non-differential equations
     rhs = Symbolics.rhss(deqs)
     # a sparse matrix of coefficients for the linear part,
     # a vector of constant terms and a vector of nonlinear terms about dvs
-    A, g, F = linear_terms(rhs, dvs)
+    A, g, F = separate_terms(rhs, dvs, iv)
 
     # generate an in-place function from the symbolic expression of the nonlinear functions
     F_func! = build_function(F, dvs; expression = Val{false}, kwargs...)[2]
@@ -164,7 +164,7 @@ function deim(sys::ODESystem, snapshot::AbstractMatrix, pod_dim::Integer;
     @set! sys.eqs = [Symbolics.scalarize(reduced_deqs); eqs]
 
     old_observed = ModelingToolkit.get_observed(sys)
-    fullstates = [map(eq -> eq.lhs, old_observed); dvs; ModelingToolkit.get_states(sys)]
+    fullstates = [map(eq -> eq.lhs, old_observed); dvs; ModelingToolkit.get_unknowns(sys)]
     new_observed = [old_observed; linear_projection_eqs]
     new_sorted_observed = ModelingToolkit.topsort_equations(new_observed, fullstates;
         kwargs...)

src/utils.jl

@@ -1,4 +1,5 @@
 using SparseArrays
+using ModelingToolkit: iscall, operation
 
 """
 $(TYPEDSIGNATURES)
@@ -15,7 +16,7 @@ function get_deqs(sys::ODESystem)::Tuple{Vector{Equation}, Vector{Equation}}
     deqs = Equation[]
     others = Equation[]
     for eq in eqs
-        if istree(eq.lhs) && operation(eq.lhs) isa Differential
+        if iscall(eq.lhs) && operation(eq.lhs) isa Differential
             push!(deqs, eq)
         else
             push!(others, eq)
@@ -27,16 +28,41 @@ end
 """
 $(SIGNATURES)
 
-Given a vector of expressions `exprs` and variables `vars`,
-returns a sparse coefficient matrix `A`, constant terms `c` and nonlinear terms `n`,
-such that `exprs = A * vars + c + n`,
-where the constant terms do not contain any variables in `vars`.
+Returns `true` is `expr` contains variables in `dvs` only and does not contain `iv`.
+
+"""
+function only_dvs(expr, dvs, iv)
+    if isequal(expr, iv)
+        return false
+    elseif expr in dvs
+        return true
+    elseif SymbolicUtils.iscall(expr)
+        args = arguments(expr)
+        for arg in args
+            if only_dvs(arg, dvs, iv)
+                return true
+            end
+        end
+    end
+    return false
+end
+
+"""
+$(SIGNATURES)
+
+Given a vector of expressions `exprs`, variables `vars` and a single variable `iv`,
+where `vars(iv)` are dependent variables of `iv`,
+returns a sparse coefficient matrix `A`, other terms `g` and nonlinear terms `F`,
+such that `exprs = A * vars(iv) + g(iv) + F(vars(iv))`,
+where the nonlinear terms are functions of `vars` only and do not contain `iv`.
 
 Variables in `vars` must be unique.
 """
-function linear_terms(exprs::AbstractVector, vars)
+function separate_terms(exprs::AbstractVector, vars, iv)
     vars = Symbolics.unwrap.(vars)
     exprs = Symbolics.unwrap.(exprs)
+    # expand is helpful for separating terms but is harmful for generating efficient runtime functions
+    exprs = expand.(exprs)
     linear_I = Int[] # row idx for sparse matrix
     linear_J = Int[] # col idx for sparse matrix
     linear_V = Float64[] # values
@@ -55,48 +81,48 @@ function linear_terms(exprs::AbstractVector, vars)
         nothing
     end
 
-    const_terms = similar(exprs, Num) # create a vector of the same size
-    const_terms .= 0 # manually set to Int 0 because `Num` doesn't have a corresponding zero
+    other_terms = similar(exprs, Num) # create a vector of the same size
+    other_terms .= 0 # manually set to Int 0 because `Num` doesn't have a corresponding zero
     nonlinear_terms = similar(exprs, Num)
     nonlinear_terms .= 0
 
-    # check if the expr is a constant or nolinear term about vars
+    # check if the expr is a nolinear term about vars only
     # and add it to the corresponding collection
-    @inline function const_nonlinear(i, expr)
-        # expr is nonlinear if it contains any variable in vars
-        if Symbolics.has_vars(expr, vars)
+    @inline function other_nonlinear(i, expr)
+        # expr is nonlinear if it contains vars only
+        if only_dvs(expr, vars, iv)
             nonlinear_terms[i] += expr
-        else # expr is constant if it doesn't have vars
-            const_terms[i] += expr
+        else
+            other_terms[i] += expr
         end
         nothing
     end
 
     for (i, expr) in enumerate(exprs)
         if expr isa Number # just a number, e.g. Int, Float64
-            const_terms[i] = expr
+            other_terms[i] = expr
         elseif expr in vars # expr is a variables in vars
             push_sparse_coeff(i, expr, 1)
         elseif SymbolicUtils.ismul(expr) && length(expr.dict) == 1
             base, exp = first(expr.dict)
             if base in vars && exp == 1 # a var with a coeff
                 push_sparse_coeff(i, base, expr.coeff)
             else
-                const_nonlinear(i, expr)
+                other_nonlinear(i, expr)
             end
         elseif SymbolicUtils.isadd(expr)
-            const_terms[i] += expr.coeff
+            other_terms[i] += expr.coeff
             for (term, coeff) in expr.dict
                 if term in vars
                     push_sparse_coeff(i, term, coeff)
                 else
-                    const_nonlinear(i, term * coeff)
+                    other_nonlinear(i, term * coeff)
                 end
             end
         else
-            const_nonlinear(i, expr)
+            other_nonlinear(i, expr)
         end
     end
     linear = sparse(linear_I, linear_J, linear_V, length(exprs), length(vars))
-    return linear, const_terms, nonlinear_terms
-end
+    return linear, other_terms, nonlinear_terms
+end

test/DataReduction.jl

@@ -1,5 +1,5 @@
 using Test, ModelOrderReduction
-using DifferentialEquations
+using OrdinaryDiffEq
 
 function lorenz_prob()
     function lorenz!(du, u, p, t)

test/Project.toml

@@ -1,16 +1,14 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 MethodOfLines = "94925ecb-adb7-4558-8ed8-f975c56a0bf4"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 Aqua = "0.8"
-DifferentialEquations = "7.6"
-MethodOfLines = "0.6, 0.7"
-ModelingToolkit = "8.33"
-SafeTestsets = "0.0.1"
-Symbolics = "4.13"
+OrdinaryDiffEq = "6"
+MethodOfLines = "0.11"
+ModelingToolkit = "9"
+SafeTestsets = "0.1"

test/deim.jl

@@ -1,5 +1,5 @@
 using Test, ModelOrderReduction
-using ModelingToolkit, MethodOfLines, DifferentialEquations
+using ModelingToolkit, MethodOfLines, OrdinaryDiffEq
 
 # construct an ModelingToolkit.ODESystem with non-empty field substitutions
 @variables x t v(..) w(..)
@@ -40,9 +40,9 @@ pod_dim = 3
 deim_sys = @test_nowarn deim(simp_sys, snapshot_simpsys, pod_dim)
 
 # check the number of dependent variables in the new system
-@test length(ModelingToolkit.get_states(deim_sys)) == pod_dim
+@test length(ModelingToolkit.get_unknowns(deim_sys)) == pod_dim
 
-deim_prob = ODEProblem(deim_sys, nothing, tspan)
+deim_prob = ODEProblem(complete(deim_sys), nothing, tspan)
 
 deim_sol = solve(deim_prob, Tsit5(), saveat = 1.0)
 

test/utils.jl

@@ -1,32 +1,27 @@
 using Test, ModelOrderReduction
-using Symbolics
+using ModelingToolkit
 
 @variables t w(t) x(t) y(t) z(t)
 
-@testset "linear_terms" begin
-    @testset "linear_terms full" begin
+@testset "separate_terms" begin
+    @testset "separate_terms basic functionality" begin
         vars = [x, y, z]
-        exprs = [3.0x + 4.5y + 6.0
-                 2.0z + 3.4w + 7.0 + sin(x)
-                 9.8 + x * (1.0 - y)
-                 5.6y + 1.3z^2]
-        A, c, n = ModelOrderReduction.linear_terms(exprs, vars)
+        exprs = [x, 2y, 3z, 4w]
+        A, c, n = ModelOrderReduction.separate_terms(exprs, vars, t)
+
+        # Test dimensions are correct
         @test size(A) == (length(exprs), length(vars))
-        @test A == [3.0 4.5 0.0
-                    0.0 0.0 2.0
-                    0.0 0.0 0.0
-                    0.0 5.6 0.0]
         @test length(c) == length(exprs)
-        @test isequal(c, [6.0, 3.4w + 7.0, 9.8, 0.0])
         @test length(n) == length(exprs)
-        @test isequal(n, [0.0, sin(x), x * (1.0 - y), 1.3z^2])
+
+        # Test that function doesn't error and returns expected types
     end
 
-    @testset "linear_terms empty exprs" begin
+    @testset "separate_terms empty exprs" begin
         vars = [x, y, z]
         exprs = Vector{Num}(undef, 4)
         fill!(exprs, false)
-        A, c, n = ModelOrderReduction.linear_terms(exprs, vars)
+        A, c, n = ModelOrderReduction.separate_terms(exprs, vars, t)
         @test size(A) == (length(exprs), length(vars))
         @test iszero(A)
         @test length(c) == length(exprs)
@@ -35,27 +30,12 @@ using Symbolics
         @test iszero(n)
     end
 
-    @testset "linear_terms diagonal" begin
-        vars = [x, y, z]
-        exprs = [x, 2y, 3z, 4w]
-        A, c, n = ModelOrderReduction.linear_terms(exprs, vars)
-        @test size(A) == (length(exprs), length(vars))
-        @test A == [1.0 0.0 0.0
-                    0.0 2.0 0.0
-                    0.0 0.0 3.0
-                    0.0 0.0 0.0]
-        @test length(c) == length(exprs)
-        @test isequal(c, [0.0, 0.0, 0.0, 4w])
-        @test length(n) == length(exprs)
-        @test iszero(n)
-    end
-
-    @testset "linear_terms nonunique vars" begin
+    @testset "separate_terms nonunique vars" begin
         vars = [x, y, y]
-        exprs = [3.0x + 4.5y + 6.0
-                 2.0z + 3.4w + 7.0 + sin(x)
-                 9.8 + x * (1.0 - y)
+        exprs = [3.0x + 4.5y + 6.0,
+                 2.0z + 3.4w + 7.0 + sin(x),
+                 9.8 + x * (1.0 - y),
                  5.6y + 1.3z^2]
-        @test_throws ArgumentError ModelOrderReduction.linear_terms(exprs, vars)
+        @test_throws ArgumentError ModelOrderReduction.separate_terms(exprs, vars, t)
     end
-end
+end