Merge branch 'master' into myb/ss

Author: YingboMa

.github/workflows/Downstream.yml

@@ -18,6 +18,7 @@ jobs:
         os: [ubuntu-latest]
         package:
           - {user: SciML, repo: Catalyst.jl, group: All}
+          - {user: SciML, repo: CellMLToolkit.jl, group: All}
           - {user: SciML, repo: NeuralPDE.jl, group: NNPDE}
           - {user: SciML, repo: DataDrivenDiffEq.jl, group: Standard}
 

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "5.13.6"
+version = "5.14.2"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -51,7 +51,7 @@ Distributions = "0.23, 0.24"
 DocStringExtensions = "0.7, 0.8"
 IfElse = "0.1"
 LabelledArrays = "1.3"
-Latexify = "0.11, 0.12, 0.13, 0.14"
+Latexify = "0.11, 0.12, 0.13, 0.14, 0.15"
 LightGraphs = "1.3"
 MacroTools = "0.5"
 NaNMath = "0.3"
@@ -65,8 +65,8 @@ SciMLBase = "1.3"
 Setfield = "0.7"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
-SymbolicUtils = "0.8.3, 0.9"
-Symbolics = "0.1.8"
+SymbolicUtils = "0.10.1"
+Symbolics = "0.1.14"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"
 julia = "1.2"
@@ -81,7 +81,8 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
+Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "ForwardDiff", "GalacticOptim", "NonlinearSolve", "OrdinaryDiffEq", "Optim", "Random", "SteadyStateDiffEq", "Test", "StochasticDiffEq"]
+test = ["BenchmarkTools", "ForwardDiff", "GalacticOptim", "NonlinearSolve", "OrdinaryDiffEq", "Optim", "Random", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]

docs/src/basics/ContextualVariables.md

@@ -7,7 +7,7 @@ the `@variable` which is defined by
 [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl). For example:
 
 ```julia
-@variabes x y(x)
+@variables x y(x)
 ```
 
 This is used for the "normal" variable of a given system, like the states of a

docs/src/tutorials/ode_modeling.md

@@ -13,7 +13,7 @@ But if you want to just see some code and run, here's an example:
 ```julia
 using ModelingToolkit
 
-@variables t, x(t) RHS(t)  # independent and dependent variables
+@variables t x(t) RHS(t)  # independent and dependent variables
 @parameters τ       # parameters
 D = Differential(t) # define an operator for the differentiation w.r.t. time
 
@@ -51,7 +51,7 @@ first set the forcing function to a constant value.
 ```julia
 using ModelingToolkit
 
-@variables t, x(t)  # independent and dependent variables
+@variables t x(t)  # independent and dependent variables
 @parameters τ       # parameters
 D = Differential(t) # define an operator for the differentiation w.r.t. time
 

docs/src/tutorials/tearing_parallelism.md

@@ -37,7 +37,7 @@ function connect_heat(ps...)
 end
 
 # Basic electric components
-const t = Sym{ModelingToolkit.Parameter{Real}}(:t)
+@parameters t
 const D = Differential(t)
 function Pin(;name)
     @variables v(t) i(t)

examples/electrical_components.jl

@@ -2,7 +2,6 @@ using Test
 using ModelingToolkit, OrdinaryDiffEq
 
 # Basic electric components
-#const t = Sym{ModelingToolkit.Parameter{Real}}(:t)
 @parameters t
 function Pin(;name)
     @variables v(t) i(t)

src/ModelingToolkit.jl

@@ -137,6 +137,7 @@ include("structural_transformation/StructuralTransformations.jl")
 @reexport using .StructuralTransformations
 
 export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr
+export DAEFunctionExpr, DAEProblemExpr
 export SDESystem, SDEFunction, SDEFunctionExpr, SDESystemExpr
 export SystemStructure
 export JumpSystem
@@ -157,6 +158,7 @@ export Differential, expand_derivatives, @derivatives
 export IntervalDomain, ProductDomain, ⊗, CircleDomain
 export Equation, ConstrainedEquation
 export Term, Sym
+export SymScope, LocalScope, ParentScope, GlobalScope
 export independent_variable, states, parameters, equations, controls, observed, structure
 export structural_simplify
 

src/parameters.jl

@@ -1,24 +1,35 @@
-import SymbolicUtils: symtype, term
-struct Parameter{T} end
+import SymbolicUtils: symtype, term, hasmetadata
+struct MTKParameterCtx end
 
+isparameter(x::Num) = isparameter(value(x))
+isparameter(x::Symbolic) = getmetadata(x, MTKParameterCtx, false)
 isparameter(x) = false
-isparameter(::Sym{<:Parameter}) = true
-isparameter(::Sym{<:FnType{<:Any, <:Parameter}}) = true
 
-SymbolicUtils.@number_methods(Sym{Parameter{Real}},
-                              term(f, a),
-                              term(f, a, b), skipbasics)
+"""
+    toparam(s::Sym)
 
-SymbolicUtils.symtype(s::Symbolic{Parameter{T}}) where T = T
-SymbolicUtils.similarterm(t::Term{<:Parameter}, f, args) = Term(f, args)
+Maps the variable to a paramter.
+"""
+toparam(s::Symbolic) = setmetadata(s, MTKParameterCtx, true)
+toparam(s::Num) = Num(toparam(value(s)))
 
-Base.convert(::Type{Num}, x::Symbolic{Parameter{T}}) where {T<:Number} = Num(x)
+"""
+    tovar(s::Sym)
+
+Maps the variable to a state.
+"""
+tovar(s::Symbolic) = setmetadata(s, MTKParameterCtx, false)
+tovar(s::Num) = Num(tovar(value(s)))
 
 """
 $(SIGNATURES)
 
 Define one or more known variables.
 """
 macro parameters(xs...)
-    esc(_parse_vars(:parameters, Parameter{Real}, xs))
+    Symbolics._parse_vars(:parameters,
+                          Real,
+                          xs,
+                          x -> x isa Array ? toparam.(x) : toparam(x)
+                         ) |> esc
 end

src/structural_transformation/codegen.jl

@@ -124,7 +124,7 @@ function partitions_dag(s::SystemStructure)
     sparse(I, J, true, n, n)
 end
 
-function gen_nlsolve(sys, eqs, vars)
+function gen_nlsolve(sys, eqs, vars; checkbounds=true)
     @assert !isempty(vars)
     @assert length(eqs) == length(vars)
     rhss = map(x->x.rhs, eqs)
@@ -142,8 +142,8 @@ function gen_nlsolve(sys, eqs, vars)
     fname = gensym("fun")
     f = Func(
         [
-         DestructuredArgs(vars)
-         DestructuredArgs(params)
+         DestructuredArgs(vars, inbounds=!checkbounds)
+         DestructuredArgs(params, inbounds=!checkbounds)
         ],
         [],
         isscalar ? rhss[1] : MakeArray(rhss, SVector)
@@ -161,7 +161,7 @@ function gen_nlsolve(sys, eqs, vars)
 
     [
      fname ← @RuntimeGeneratedFunction(f)
-     DestructuredArgs(vars) ← solver_call
+     DestructuredArgs(vars, inbounds=!checkbounds) ← solver_call
     ]
 end
 
@@ -205,8 +205,8 @@ function build_torn_function(
         Func(
              [
               out
-              DestructuredArgs(states)
-              DestructuredArgs(parameters(sys))
+              DestructuredArgs(states, inbounds=!checkbounds)
+              DestructuredArgs(parameters(sys), inbounds=!checkbounds)
               independent_variable(sys)
              ],
              [],
@@ -222,7 +222,7 @@ function build_torn_function(
         observedfun = let sys = sys, dict = Dict()
             function generated_observed(obsvar, u, p, t)
                 obs = get!(dict, value(obsvar)) do
-                    build_observed_function(sys, obsvar)
+                    build_observed_function(sys, obsvar, checkbounds=checkbounds)
                 end
                 obs(u, p, t)
             end
@@ -257,7 +257,8 @@ end
 function build_observed_function(
         sys, syms;
         expression=false,
-        output_type=Array
+        output_type=Array,
+        checkbounds=true
     )
 
     if (isscalar = !(syms isa Vector))
@@ -293,7 +294,7 @@ function build_observed_function(
         torn_eqs  = map(idxs-> eqs[idxs.e_residual], subset)
         torn_vars = map(idxs->fullvars[idxs.v_residual], subset)
 
-        solves = gen_nlsolve.((sys,), torn_eqs, torn_vars)
+        solves = gen_nlsolve.((sys,), torn_eqs, torn_vars; checkbounds=checkbounds)
     else
         solves = []
     end
@@ -308,8 +309,8 @@ function build_observed_function(
 
     ex = Func(
         [
-         DestructuredArgs(diffvars)
-         DestructuredArgs(parameters(sys))
+         DestructuredArgs(diffvars, inbounds=!checkbounds)
+         DestructuredArgs(parameters(sys), inbounds=!checkbounds)
          independent_variable(sys)
         ],
         [],

src/structural_transformation/pantelides.jl

@@ -2,23 +2,6 @@
 ### Reassemble: structural information -> system
 ###
 
-# Naive subtree matching, we can make the dict have levels
-# Going forward we should look into storing the depth in Terms
-function walk_and_substitute(expr, substitution_dict)
-    if haskey(substitution_dict, expr)
-        return substitution_dict[expr]
-    elseif istree(expr)
-        return similarterm(
-                           expr,
-                           walk_and_substitute(operation(expr), substitution_dict),
-                           map(x->walk_and_substitute(x, substitution_dict), arguments(expr))
-                          )
-    else
-        @assert !(expr isa Equation) "RHS cannot contain equations"
-        return expr
-    end
-end
-
 function pantelides_reassemble(sys, eqassoc, assign)
     s = structure(sys)
     @unpack fullvars, varassoc = s
@@ -74,7 +57,7 @@ function pantelides_reassemble(sys, eqassoc, assign)
         end
         rhs = ModelingToolkit.expand_derivatives(D(eq.rhs))
         substitution_dict = Dict(x.lhs => x.rhs for x in out_eqs if x !== nothing && x.lhs isa Symbolic)
-        sub_rhs = walk_and_substitute(rhs, substitution_dict)
+        sub_rhs = substitute(rhs, substitution_dict)
         out_eqs[e] = lhs ~ sub_rhs
     end