Update dependencies versions and fix CI

.github/workflows/Tests.yml

@@ -24,7 +24,6 @@ jobs:
  matrix:
  version:
  - "1"
- - "1.6"
  uses: "SciML/.github/.github/workflows/tests.yml@v1"
  with:
  julia-version: "${{ matrix.version }}"

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/make.jl

@@ -6,12 +6,12 @@ cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true)
 include("pages.jl")
 
 makedocs(sitename = "ModelOrderReduction.jl",
-  authors = "Bowen S. Zhu",
-  modules = [ModelOrderReduction],
-  clean = true, doctest = false, linkcheck = true,
-  warnonly = [:missing_docs, :example_block],
-  format = Documenter.HTML(assets = ["assets/favicon.ico"],
-  canonical = "https://docs.sciml.ai/ModelOrderReduction/stable/"),
-  pages = pages)
+ authors = "Bowen S. Zhu",
+ modules = [ModelOrderReduction],
+ clean = true, doctest = false, linkcheck = true,
+ warnonly = [:missing_docs, :example_block],
+ format = Documenter.HTML(assets = ["assets/favicon.ico"],
+ canonical = "https://docs.sciml.ai/ModelOrderReduction/stable/"),
+ pages = pages)
 
 deploydocs(repo = "github.com/SciML/ModelOrderReduction.jl.git")

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"]
 ]

docs/src/tutorials/deim_FitzHugh-Nagumo.md

@@ -107,7 +107,7 @@ svdval_w = svdvals(snapshot_w)
 svdval_fv = svdvals(snapshot_fv)
 using Plots, LaTeXStrings
 svd_plt = plot(yscale = :log10, xticks = eachindex(svdval_v), titlefont = 11,
-  legendfont = 10, title = "Singular values of the snapshots")
+ legendfont = 10, title = "Singular values of the snapshots")
 plot!(svd_plt, svdval_v, markershape = :circle, label = L"Singular Val of $v$")
 plot!(svd_plt, svdval_w, markershape = :circle, label = L"Singular Val of $w$")
 plot!(svd_plt, svdval_fv, markershape = :circle, label = L"Singular Val of $f(v)$")
@@ -118,9 +118,9 @@ points ``x`` from the full-order system.
 
 ```@example deim_FitzHugh_Nagumo
 full_plt = plot(xlabel = L"v(x,t)", ylabel = L"x", zlabel = L"w(x,t)", xlims = (-0.5, 2.0),
-  ylims = (0.0, L), zlims = (0.0, 0.25), legend = false, xflip = true,
-  camera = (50, 30), titlefont = 10,
-  title = "Phase−Space diagram of full $(nameof(pde_sys)) system")
+ ylims = (0.0, L), zlims = (0.0, 0.25), legend = false, xflip = true,
+ camera = (50, 30), titlefont = 10,
+ title = "Phase−Space diagram of full $(nameof(pde_sys)) system")
 @views for i in 1:nₓ
  plot!(full_plt, snapshot_v[i, :], _ -> sol_x[i], snapshot_w[i, :])
 end
@@ -149,9 +149,9 @@ And plot the result from the POD-DEIM reduced system.
 
 ```@example deim_FitzHugh_Nagumo
 deim_plt = plot(xlabel = L"v(x,t)", ylabel = L"x", zlabel = L"w(x,t)", xlims = (-0.5, 2.0),
-  ylims = (0.0, L), zlims = (0.0, 0.25), legend = false, xflip = true,
-  camera = (50, 30), titlefont = 10,
-  title = "Phase−Space diagram of reduced $(nameof(pde_sys)) system")
+ ylims = (0.0, L), zlims = (0.0, 0.25), legend = false, xflip = true,
+ camera = (50, 30), titlefont = 10,
+ title = "Phase−Space diagram of reduced $(nameof(pde_sys)) system")
 @views for i in 1:nₓ
  plot!(deim_plt, sol_deim_v[i, :], _ -> sol_deim_x[i], sol_deim_w[i, :])
 end
@@ -178,12 +178,12 @@ function unconnected(v::AbstractVector, nₜ::Integer)
  vec(data)
 end
 plt_2 = plot(xlabel = L"v(x,t)", ylabel = L"x", zlabel = L"w(x,t)", xlims = (-0.5, 2.0),
-  ylims = (0.0, L), zlims = (0.0, 0.25), xflip = true, camera = (50, 30),
-  titlefont = 10, title = "Comparison of full and reduced systems")
+ ylims = (0.0, L), zlims = (0.0, 0.25), xflip = true, camera = (50, 30),
+ titlefont = 10, title = "Comparison of full and reduced systems")
 plot!(plt_2, unconnected(snapshot_v), unconnected(sol_x, nₜ), unconnected(snapshot_w),
-  label = "Full$(length(ode_sys.eqs))")
+ label = "Full$(length(ode_sys.eqs))")
 plot!(plt_2, unconnected(sol_deim_v), unconnected(sol_deim_x, nₜ_deim),
-  unconnected(sol_deim_w), label = "POD$(pod_dim)/DEIM$(deim_dim)")
+ unconnected(sol_deim_w), label = "POD$(pod_dim)/DEIM$(deim_dim)")
 ```
 
 As we can see, the reduced-order system captures the limit cycle of the original full-order

src/DataReduction/POD.jl

@@ -39,9 +39,9 @@ mutable struct POD <: AbstractDRProblem
  spectrum::Any
  # constructors
  function POD(snaps;
-  min_renergy = 1.0,
-  min_nmodes::Int = 1,
-  max_nmodes::Int = length(snaps[1]))
+ min_renergy = 1.0,
+ min_nmodes::Int = 1,
+ max_nmodes::Int = length(snaps[1]))
  nmodes = min_nmodes
  errorhandle(snaps, nmodes, min_renergy, min_nmodes, max_nmodes)
  new(snaps, min_renergy, min_nmodes, max_nmodes, nmodes, missing, 1.0, missing)
@@ -66,7 +66,7 @@ end
 function reduce!(pod::POD, alg::SVD)
  u, s, v = _svd(pod.snapshots; alg.kwargs...)
  pod.nmodes, pod.renergy = determine_truncation(s, pod.min_nmodes, pod.max_nmodes,
-  pod.min_renergy)
+ pod.min_renergy)
  pod.rbasis = u[:, 1:(pod.nmodes)]
  pod.spectrum = s
  nothing
@@ -94,6 +94,6 @@ function Base.show(io::IO, pod::POD)
  print(io, "POD \n")
  print(io, "Reduction Order = ", pod.nmodes, "\n")
  print(io, "Snapshot size = (", size(pod.snapshots, 1), ",", size(pod.snapshots[1], 2),
-  ")\n")
+ ")\n")
  print(io, "Relative Energy = ", pod.renergy, "\n")
 end

src/ErrorHandle.jl

@@ -7,7 +7,7 @@ function errorhandle(data::AbstractMatrix, nmodes::Int, max_energy, min_nmodes,
 end
 
 function errorhandle(data::AbstractVector{T}, nmodes::Int, max_energy, min_nmodes,
-  max_nmodes) where {T <: AbstractVector}
+ max_nmodes) where {T <: AbstractVector}
  @assert size(data[1], 1)>1 "State vector is expected to be vector valued."
  s = min(size(data, 1), size(data[1], 1))
  @assert 0<nmodes<=s "Number of modes should be in {1,2,...,$s}."

src/ModelOrderReduction.jl

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

src/deim.jl

@@ -65,9 +65,9 @@ the ``\\rho_i``-th column of the identity matrix ``I_n\\in\\mathbb R^{n\\times n
 - `linear_projection_eqs`: the linear projection mapping ``\\mathbf y=V\\hat{\\mathbf y}``.
 """
 function deim(full_vars::AbstractVector, linear_coeffs::AbstractMatrix,
-  constant_part::AbstractVector, nonlinear_part::AbstractVector,
-  reduced_vars::AbstractVector, linear_projection_matrix::AbstractMatrix,
-  nonlinear_projection_matrix::AbstractMatrix; kwargs...)
+ constant_part::AbstractVector, nonlinear_part::AbstractVector,
+ reduced_vars::AbstractVector, linear_projection_matrix::AbstractMatrix,
+ nonlinear_projection_matrix::AbstractMatrix; kwargs...)
  # rename variables for convenience
  y = full_vars
  A = linear_coeffs
@@ -119,8 +119,8 @@ nonlinear terms, which is computed by executing the runtime-generated function f
 nonlinear expressions.
 """
 function deim(sys::ODESystem, snapshot::AbstractMatrix, pod_dim::Integer;
-  deim_dim::Integer = pod_dim, name::Symbol = Symbol(nameof(sys), :_deim),
-  kwargs...)::ODESystem
+ deim_dim::Integer = pod_dim, name::Symbol = Symbol(nameof(sys), :_deim),
+ kwargs...)::ODESystem
  sys = deepcopy(sys)
  @set! sys.name = name
 
@@ -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,10 +164,10 @@ 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...)
+ kwargs...)
  @set! sys.observed = new_sorted_observed
 
  inv_dict = Dict(Symbolics.scalarize(ŷ .=> V' * dvs)) # reduced vars to original vars

src/utils.jl

@@ -27,16 +27,41 @@
 """
 $(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.istree(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 +80,48 @@
  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
+ return linear, other_terms, nonlinear_terms
 end

test/Project.toml

@@ -4,13 +4,11 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 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"
+DifferentialEquations = "7.13"
+MethodOfLines = "0.11"
+ModelingToolkit = "9"
+SafeTestsets = "0.1"

test/deim.jl

@@ -40,7 +40,7 @@ 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)
 

test/runtests.jl

@@ -1,6 +1,14 @@
 using SafeTestsets
 
-@safetestset "Quality Assurance" begin include("qa.jl") end
-@safetestset "POD" begin include("DataReduction.jl") end
-@safetestset "utils" begin include("utils.jl") end
-@safetestset "DEIM" begin include("deim.jl") end
+@safetestset "Quality Assurance" begin
+ include("qa.jl")
+end
+@safetestset "POD" begin
+ include("DataReduction.jl")
+end
+@safetestset "utils" begin
+ include("utils.jl")
+end
+@safetestset "DEIM" begin
+ include("deim.jl")
+end