Merge origin/master

Author: chriselrod

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/typos@v1.16.25 
+        uses: crate-ci/typos@v1.17.0 

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <[email protected]>", "Chris Rackauckas <[email protected]> and contributors"]
-version = "8.74.1"
+version = "9.0.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -17,6 +17,7 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
+DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 FunctionWrappersWrappers = "77dc65aa-8811-40c2-897b-53d922fa7daf"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
@@ -72,6 +73,7 @@ Distributed = "1"
 Distributions = "0.23, 0.24, 0.25"
 DocStringExtensions = "0.7, 0.8, 0.9"
 DomainSets = "0.6"
+DynamicQuantities = "0.8, 0.9, 0.10"
 ForwardDiff = "0.10.3"
 FunctionWrappersWrappers = "0.1"
 Graphs = "1.5.2"
@@ -98,7 +100,7 @@ SimpleNonlinearSolve = "0.1.0, 1"
 SparseArrays = "1"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
-SymbolicIndexingInterface = "0.3"
+SymbolicIndexingInterface = "0.3.1"
 SymbolicUtils = "1.0"
 Symbolics = "5.7"
 URIs = "1"

docs/Project.toml

@@ -27,14 +27,14 @@ Distributions = "0.25"
 Documenter = "1"
 ModelingToolkit = "8.33"
 ModelingToolkitDesigner = "1"
-NonlinearSolve = "0.3, 1, 2"
+NonlinearSolve = "0.3, 1, 2, 3"
 Optim = "1.7"
 Optimization = "3.9"
 OptimizationOptimJL = "0.1"
 OrdinaryDiffEq = "6.31"
 Plots = "1.36"
 StochasticDiffEq = "6"
-StructuralIdentifiability = "0.4"
+StructuralIdentifiability = "0.4, 0.5"
 SymbolicUtils = "1"
 Symbolics = "5"
 Unitful = "1.12"

ext/MTKBifurcationKitExt.jl

@@ -34,11 +34,11 @@ struct ObservableRecordFromSolution{S, T}
         param_end_idxs = state_end_idxs + length(parameters(nsys))
 
         bif_par_idx = state_end_idxs + bif_idx
-        # Gets the (base) substitution values for states. 
+        # Gets the (base) substitution values for states.
         subs_vals_states = Pair.(states(nsys), u0_vals)
-        # Gets the (base) substitution values for parameters. 
+        # Gets the (base) substitution values for parameters.
         subs_vals_params = Pair.(parameters(nsys), p_vals)
-        # Gets the (base) substitution values for observables. 
+        # Gets the (base) substitution values for observables.
         subs_vals_obs = [obs.lhs => substitute(obs.rhs,
             [subs_vals_states; subs_vals_params]) for obs in observed(nsys)]
         # Sometimes observables depend on other observables, hence we make a second update to this vector.
@@ -136,7 +136,7 @@ function BifurcationKit.BifurcationProblem(osys::ODESystem, args...; kwargs...)
     nsys = NonlinearSystem([0 ~ eq.rhs for eq in equations(osys)],
         states(osys),
         parameters(osys);
-        name = osys.name)
+        name = nameof(osys))
     return BifurcationKit.BifurcationProblem(nsys, args...; kwargs...)
 end
 

format/Project.toml

@@ -0,0 +1,2 @@
+[deps]
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"

src/ModelingToolkit.jl

@@ -152,6 +152,7 @@ include("systems/pde/pdesystem.jl")
 
 include("systems/sparsematrixclil.jl")
 include("systems/discrete_system/discrete_system.jl")
+include("systems/unit_check.jl")
 include("systems/validation.jl")
 include("systems/dependency_graphs.jl")
 include("clock.jl")

src/parameters.jl

@@ -100,6 +100,11 @@ function split_parameters_by_type(ps)
         end
         tighten_types = x -> identity.(x)
         split_ps = tighten_types.(Base.Fix1(getindex, ps).(split_idxs))
+
+        if ps isa StaticArray
+            parrs = map(x -> SArray{Tuple{size(x)...}}(x), split_ps)
+            split_ps = SArray{Tuple{size(parrs)...}}(parrs)
+        end
         if length(split_ps) == 1  #Tuple not needed, only 1 type
             return split_ps[1], split_idxs
         else

src/systems/abstractsystem.jl

@@ -185,75 +185,81 @@ end
 #Treat the result as a vector of symbols always
 function SymbolicIndexingInterface.is_variable(sys::AbstractSystem, sym)
     if unwrap(sym) isa Int    # [x, 1] coerces 1 to a Num
-        return unwrap(sym) in 1:length(unknown_states(sys))
+        return unwrap(sym) in 1:length(variable_symbols(sys))
     end
-    return any(isequal(sym), unknown_states(sys)) ||
+    return any(isequal(sym), variable_symbols(sys)) ||
            hasname(sym) && is_variable(sys, getname(sym))
 end
 
 function SymbolicIndexingInterface.is_variable(sys::AbstractSystem, sym::Symbol)
-    return any(isequal(sym), getname.(unknown_states(sys))) ||
+    return any(isequal(sym), getname.(variable_symbols(sys))) ||
            count('₊', string(sym)) == 1 &&
-           count(isequal(sym), Symbol.(sys.name, :₊, getname.(unknown_states(sys)))) == 1
+           count(isequal(sym), Symbol.(nameof(sys), :₊, getname.(variable_symbols(sys)))) ==
+           1
 end
 
 function SymbolicIndexingInterface.variable_index(sys::AbstractSystem, sym)
     if unwrap(sym) isa Int
         return unwrap(sym)
     end
-    idx = findfirst(isequal(sym), unknown_states(sys))
+    idx = findfirst(isequal(sym), variable_symbols(sys))
     if idx === nothing && hasname(sym)
         idx = variable_index(sys, getname(sym))
     end
     return idx
 end
 
 function SymbolicIndexingInterface.variable_index(sys::AbstractSystem, sym::Symbol)
-    idx = findfirst(isequal(sym), getname.(unknown_states(sys)))
+    idx = findfirst(isequal(sym), getname.(variable_symbols(sys)))
     if idx !== nothing
         return idx
     elseif count('₊', string(sym)) == 1
-        return findfirst(isequal(sym), Symbol.(sys.name, :₊, getname.(unknown_states(sys))))
+        return findfirst(isequal(sym),
+            Symbol.(nameof(sys), :₊, getname.(variable_symbols(sys))))
     end
     return nothing
 end
 
+SymbolicIndexingInterface.variable_symbols(sys::AbstractMultivariateSystem) = sys.dvs
+
 function SymbolicIndexingInterface.variable_symbols(sys::AbstractSystem)
     return unknown_states(sys)
 end
 
 function SymbolicIndexingInterface.is_parameter(sys::AbstractSystem, sym)
     if unwrap(sym) isa Int
-        return unwrap(sym) in 1:length(parameters(sys))
+        return unwrap(sym) in 1:length(parameter_symbols(sys))
     end
 
-    return any(isequal(sym), parameters(sys)) ||
+    return any(isequal(sym), parameter_symbols(sys)) ||
            hasname(sym) && is_parameter(sys, getname(sym))
 end
 
 function SymbolicIndexingInterface.is_parameter(sys::AbstractSystem, sym::Symbol)
-    return any(isequal(sym), getname.(parameters(sys))) ||
+    return any(isequal(sym), getname.(parameter_symbols(sys))) ||
            count('₊', string(sym)) == 1 &&
-           count(isequal(sym), Symbol.(sys.name, :₊, getname.(parameters(sys)))) == 1
+           count(isequal(sym),
+        Symbol.(nameof(sys), :₊, getname.(parameter_symbols(sys)))) == 1
 end
 
 function SymbolicIndexingInterface.parameter_index(sys::AbstractSystem, sym)
     if unwrap(sym) isa Int
         return unwrap(sym)
     end
-    idx = findfirst(isequal(sym), parameters(sys))
+    idx = findfirst(isequal(sym), parameter_symbols(sys))
     if idx === nothing && hasname(sym)
         idx = parameter_index(sys, getname(sym))
     end
     return idx
 end
 
 function SymbolicIndexingInterface.parameter_index(sys::AbstractSystem, sym::Symbol)
-    idx = findfirst(isequal(sym), getname.(parameters(sys)))
+    idx = findfirst(isequal(sym), getname.(parameter_symbols(sys)))
     if idx !== nothing
         return idx
     elseif count('₊', string(sym)) == 1
-        return findfirst(isequal(sym), Symbol.(sys.name, :₊, getname.(parameters(sys))))
+        return findfirst(isequal(sym),
+            Symbol.(nameof(sys), :₊, getname.(parameter_symbols(sys))))
     end
     return nothing
 end
@@ -263,7 +269,7 @@ function SymbolicIndexingInterface.parameter_symbols(sys::AbstractSystem)
 end
 
 function SymbolicIndexingInterface.is_independent_variable(sys::AbstractSystem, sym)
-    return any(isequal(sym), independent_variables(sys))
+    return any(isequal(sym), independent_variable_symbols(sys))
 end
 
 function SymbolicIndexingInterface.is_independent_variable(sys::AbstractSystem, sym::Symbol)
@@ -284,6 +290,20 @@ SymbolicIndexingInterface.is_time_dependent(::AbstractTimeIndependentSystem) = f
 
 SymbolicIndexingInterface.constant_structure(::AbstractSystem) = true
 
+function SymbolicIndexingInterface.all_variable_symbols(sys::AbstractSystem)
+    syms = variable_symbols(sys)
+    obs = getproperty.(observed(sys), :lhs)
+    return isempty(obs) ? syms : vcat(syms, obs)
+end
+
+function SymbolicIndexingInterface.all_symbols(sys::AbstractSystem)
+    syms = all_variable_symbols(sys)
+    for other in (parameter_symbols(sys), independent_variable_symbols(sys))
+        isempty(other) || (syms = vcat(syms, other))
+    end
+    return syms
+end
+
 iscomplete(sys::AbstractSystem) = isdefined(sys, :complete) && getfield(sys, :complete)
 
 """
@@ -303,6 +323,7 @@ for prop in [:eqs
     :states
     :ps
     :tspan
+    :name
     :var_to_name
     :ctrls
     :defaults
@@ -639,7 +660,7 @@ function states(sys::AbstractSystem)
     end
     isempty(nonunique_states) && return nonunique_states
     # `Vector{Any}` is incompatible with the `SymbolicIndexingInterface`, which uses
-    # `elsymtype = symbolic_type(eltype(_arg))` 
+    # `elsymtype = symbolic_type(eltype(_arg))`
     # which inappropriately returns `NotSymbolic()`
     if nonunique_states isa Vector{Any}
         nonunique_states = _nonum.(nonunique_states)
@@ -650,6 +671,12 @@ end
 
 function parameters(sys::AbstractSystem)
     ps = get_ps(sys)
+    if ps == SciMLBase.NullParameters()
+        return []
+    end
+    if eltype(ps) <: Pair
+        ps = first.(ps)
+    end
     systems = get_systems(sys)
     unique(isempty(systems) ? ps : [ps; reduce(vcat, namespace_parameters.(systems))])
 end

src/systems/diffeqs/abstractodesystem.jl

@@ -912,6 +912,13 @@ function DiffEqBase.ODEProblem(sys::AbstractODESystem, args...; kwargs...)
     ODEProblem{true}(sys, args...; kwargs...)
 end
 
+function DiffEqBase.ODEProblem(sys::AbstractODESystem,
+        u0map::StaticArray,
+        args...;
+        kwargs...)
+    ODEProblem{false, SciMLBase.FullSpecialize}(sys, u0map, args...; kwargs...)
+end
+
 function DiffEqBase.ODEProblem{true}(sys::AbstractODESystem, args...; kwargs...)
     ODEProblem{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...)
 end

src/systems/diffeqs/odesystem.jl

@@ -163,7 +163,8 @@ struct ODESystem <: AbstractODESystem
             check_equations(equations(cevents), iv)
         end
         if checks == true || (checks & CheckUnits) > 0
-            all_dimensionless([dvs; ps; iv]) || check_units(deqs)
+            u = __get_unit_type(dvs, ps, iv)
+            check_units(u, deqs)
         end
         new(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac,
             ctrl_jac, Wfact, Wfact_t, name, systems, defaults, torn_matching,

src/systems/diffeqs/sdesystem.jl

@@ -134,7 +134,8 @@ struct SDESystem <: AbstractODESystem
             check_equations(equations(cevents), iv)
         end
         if checks == true || (checks & CheckUnits) > 0
-            all_dimensionless([dvs; ps; iv]) || check_units(deqs, neqs)
+            u = __get_unit_type(dvs, ps, iv)
+            check_units(u, deqs, neqs)
         end
         new(tag, deqs, neqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac,
             ctrl_jac,

src/systems/discrete_system/discrete_system.jl

@@ -103,7 +103,8 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
             check_parameters(ps, iv)
         end
         if checks == true || (checks & CheckUnits) > 0
-            all_dimensionless([dvs; ps; iv; ctrls]) || check_units(discreteEqs)
+            u = __get_unit_type(dvs, ps, iv, ctrls)
+            check_units(u, discreteEqs)
         end
         new(tag, discreteEqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, name,
             systems,

src/systems/jumps/jumpsystem.jl

@@ -112,7 +112,8 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem
             check_parameters(ps, iv)
         end
         if checks == true || (checks & CheckUnits) > 0
-            all_dimensionless([states; ps; iv]) || check_units(ap, iv)
+            u = __get_unit_type(states, ps, iv)
+            check_units(u, ap, iv)
         end
         new{U}(tag, ap, iv, states, ps, var_to_name, observed, name, systems, defaults,
             connector_type, devents, metadata, gui_metadata, complete)

src/systems/model_parsing.jl

@@ -178,15 +178,10 @@ function parse_variable_def!(dict, mod, arg, varclass, kwargs;
             (var, def)
         end
         Expr(:ref, a, b...) => begin
+            indices = map(i -> UnitRange(i.args[2], i.args[end]), b)
             parse_variable_def!(dict, mod, a, varclass, kwargs;
-                def, indices = [eval.(b)...])
+                def, indices)
         end
-        #= Expr(:if, condition, a) => begin
-        var, def = [], []
-            for var_def in a.args
-                parse_variable_def!(dict, mod, var_def, varclass, kwargs)
-            end
-        end =#
         _ => error("$arg cannot be parsed")
     end
 end
@@ -301,7 +296,7 @@ function parse_model!(exprs, comps, ext, eqs, icon, vs, ps, sps,
         parse_equations!(exprs, eqs, dict, body)
     elseif mname == Symbol("@icon")
         isassigned(icon) && error("This model has more than one icon.")
-        parse_icon!(icon, dict, body)
+        parse_icon!(body, dict, icon, mod)
     else
         error("$mname is not handled.")
     end
@@ -614,7 +609,7 @@ function parse_equations!(exprs, eqs, dict, body)
     end
 end
 
-function parse_icon!(icon, dict, body::String)
+function parse_icon!(body::String, dict, icon, mod)
     icon_dir = get(ENV, "MTK_ICONS_DIR", joinpath(DEPOT_PATH[1], "mtk_icons"))
     dict[:icon] = icon[] = if isfile(body)
         URI("file:///" * abspath(body))
@@ -633,8 +628,8 @@ function parse_icon!(icon, dict, body::String)
     end
 end
 
-function parse_icon!(icon, dict, body::Expr)
-    parse_icon!(icon, dict, eval(body))
+function parse_icon!(body::Symbol, dict, icon, mod)
+    parse_icon!(getfield(mod, body), dict, icon, mod)
 end
 
 ### Parsing Components:

src/systems/nonlinear/nonlinearsystem.jl

@@ -88,7 +88,8 @@ struct NonlinearSystem <: AbstractTimeIndependentSystem
             tearing_state = nothing, substitutions = nothing,
             complete = false, parent = nothing; checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckUnits) > 0
-            all_dimensionless([states; ps]) || check_units(eqs)
+            u = __get_unit_type(states, ps)
+            check_units(u, eqs)
         end
         new(tag, eqs, states, ps, var_to_name, observed, jac, name, systems, defaults,
             connector_type, metadata, gui_metadata, tearing_state, substitutions, complete,

src/systems/optimization/constraints_system.jl

@@ -77,7 +77,8 @@ struct ConstraintsSystem <: AbstractTimeIndependentSystem
             tearing_state = nothing, substitutions = nothing;
             checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckUnits) > 0
-            all_dimensionless([states; ps]) || check_units(constraints)
+            u = __get_unit_type(states, ps)
+            check_units(u, constraints)
         end
         new(tag, constraints, states, ps, var_to_name, observed, jac, name, systems,
             defaults,

src/systems/optimization/optimizationsystem.jl

@@ -66,9 +66,10 @@ struct OptimizationSystem <: AbstractOptimizationSystem
             gui_metadata = nothing, complete = false, parent = nothing;
             checks::Union{Bool, Int} = true)
         if checks == true || (checks & CheckUnits) > 0
-            unwrap(op) isa Symbolic && check_units(op)
-            check_units(observed)
-            all_dimensionless([states; ps]) || check_units(constraints)
+            u = __get_unit_type(states, ps)
+            unwrap(op) isa Symbolic && check_units(u, op)
+            check_units(u, observed)
+            check_units(u, constraints)
         end
         new(tag, op, states, ps, var_to_name, observed,
             constraints, name, systems, defaults, metadata, gui_metadata, complete,