Merge pull request #943 from SciML/s/param

Use metadata to denote parameters

Author: YingboMa

Project.toml

@@ -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, 0.10"
-Symbolics = "0.1.8"
+SymbolicUtils = "0.10"
+Symbolics = "0.1.14"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"
 julia = "1.2"

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)

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/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
 

src/structural_transformation/tearing.jl

@@ -50,7 +50,7 @@ function tearing_reassemble(sys; simplify=false)
             var = fullvars[iv]
             rhs = value(solve_for(eq, var; simplify=simplify, check=false))
             # if we don't simplify the rhs and the `eq` is not solved properly
-            (!simplify && var in vars(rhs)) && (rhs = SymbolicUtils.polynormalize(rhs))
+            (!simplify && occursin(rhs, var)) && (rhs = SymbolicUtils.polynormalize(rhs))
             # Since we know `eq` is linear wrt `var`, so the round off must be a
             # linear term. We can correct the round off error by a linear
             # correction.

src/systems/abstractsystem.jl

@@ -271,10 +271,8 @@ end
 function renamespace(namespace, x)
     if x isa Num
         renamespace(namespace, value(x))
-    elseif istree(x)
-        renamespace(namespace, operation(x))(arguments(x)...)
-    elseif x isa Sym
-        Sym{symtype(x)}(renamespace(namespace,nameof(x)))
+    elseif x isa Symbolic
+        rename(x, renamespace(namespace, getname(x)))
     else
         Symbol(namespace,:₊,x)
     end
@@ -311,8 +309,7 @@ function namespace_expr(O,name,iv) where {T}
     if istree(O)
         renamed = map(a->namespace_expr(a,name,iv), arguments(O))
         if operation(O) isa Sym
-            renamed_op = rename(operation(O),renamespace(name,nameof(operation(O))))
-            Term{_symparam(O)}(renamed_op,renamed)
+            renamed_op = rename(O,renamespace(name, getname(O)))
         else
             similarterm(O,operation(O),renamed)
         end

src/systems/diffeqs/odesystem.jl

@@ -106,7 +106,7 @@ iv_from_nested_derivative(x::Term) = operation(x) isa Differential ? iv_from_nes
 iv_from_nested_derivative(x::Sym) = x
 iv_from_nested_derivative(x) = missing
 
-vars(x::Sym) = [x]
+vars(x::Sym) = Set([x])
 vars(exprs::Symbolic) = vars([exprs])
 vars(exprs) = foldl(vars!, exprs; init = Set())
 vars!(vars, eq::Equation) = (vars!(vars, eq.lhs); vars!(vars, eq.rhs); vars)

src/utils.jl

@@ -72,23 +72,6 @@ function states_to_sym(states::Set)
 end
 states_to_sym(states) = states_to_sym(Set(states))
 
-"""
-    toparam(s::Sym) -> Sym{<:Parameter}
-
-Maps the variable to a paramter.
-"""
-toparam(s::Sym) = Sym{Parameter{symtype(s)}}(s.name)
-toparam(s::Sym{<:Parameter}) = s
-
-"""
-    tovar(s::Sym) -> Sym{Real}
-    tovar(s::Sym{<:Parameter}) -> Sym{Real}
-
-Maps the variable to a state.
-"""
-tovar(s::Sym{<:Parameter}) = Sym{symtype(s)}(s.name)
-tovar(s::Sym) = s
-
 function todict(d)
     eltype(d) <: Pair || throw(ArgumentError("The variable-value mapping must be a Dict."))
     d isa Dict ? d : Dict(d)

test/rc_model.jl

@@ -2,7 +2,7 @@ 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)
     ODESystem(Equation[], t, [v, i], [], name=name, defaults=[v=>1.0, i=>1.0])

test/reactionsystem_components.jl

@@ -2,14 +2,15 @@ using ModelingToolkit, LinearAlgebra, OrdinaryDiffEq, Test
 MT = ModelingToolkit
 
 # Repressilator model
-@parameters t α₀ α K n δ β μ 
+@parameters t α₀ α K n δ β μ
 @variables m(t) P(t) R(t)
-rxs = [Reaction(α₀, nothing, [m]),
+rxs = [
+       Reaction(α₀, nothing, [m]),
        Reaction(α / (1 + (R/K)^n), nothing, [m]),
        Reaction(δ, [m], nothing),
        Reaction(β, [m], [m,P]),
-       Reaction(μ, [P], nothing)    
-]
+       Reaction(μ, [P], nothing)
+      ]
 
 specs = [m,P,R]
 pars  = [α₀,α,K,n,δ,β,μ]
@@ -19,31 +20,31 @@ pars  = [α₀,α,K,n,δ,β,μ]
 @named os₁ = convert(ODESystem, rs)
 @named os₂ = convert(ODESystem, rs)
 @named os₃ = convert(ODESystem, rs)
-connections = [os₁.R ~ os₃.P, 
+connections = [os₁.R ~ os₃.P,
                os₂.R ~ os₁.P,
                os₃.R ~ os₂.P]
 @named connected = ODESystem(connections, t, [], [], systems=[os₁,os₂,os₃])
 oderepressilator = structural_simplify(connected)
 
-pvals = [os₁.α₀ => 5e-4, 
-         os₁.α => .5, 
-         os₁.K => 40.0, 
-         os₁.n => 2, 
-         os₁.δ => (log(2)/120), 
+pvals = [os₁.α₀ => 5e-4,
+         os₁.α => .5,
+         os₁.K => 40.0,
+         os₁.n => 2,
+         os₁.δ => (log(2)/120),
          os₁.β => (20*log(2)/120),
          os₁.μ => (log(2)/600),
-         os₂.α₀ => 5e-4, 
-         os₂.α => .5, 
-         os₂.K => 40.0, 
-         os₂.n => 2, 
-         os₂.δ => (log(2)/120), 
+         os₂.α₀ => 5e-4,
+         os₂.α => .5,
+         os₂.K => 40.0,
+         os₂.n => 2,
+         os₂.δ => (log(2)/120),
          os₂.β => (20*log(2)/120),
          os₂.μ => (log(2)/600),
-         os₃.α₀ => 5e-4, 
-         os₃.α => .5, 
-         os₃.K => 40.0, 
-         os₃.n => 2, 
-         os₃.δ => (log(2)/120), 
+         os₃.α₀ => 5e-4,
+         os₃.α => .5,
+         os₃.K => 40.0,
+         os₃.n => 2,
+         os₃.δ => (log(2)/120),
          os₃.β => (20*log(2)/120),
          os₃.μ => (log(2)/600)]
 u₀    = [os₁.m => 0.0, os₁.P => 20.0, os₂.m => 0.0, os₂.P => 0.0, os₃.m => 0.0, os₃.P => 0.0]
@@ -58,8 +59,8 @@ function repress!(f, y, p, t)
     f[2] = α / (1 + (y[4] / K)^n) - δ * y[2] + α₀
     f[3] = α / (1 + (y[5] / K)^n) - δ * y[3] + α₀
     f[4] = β * y[1] - μ * y[4]
-    f[5] = β * y[2] - μ * y[5]    
-    f[6] = β * y[3] - μ * y[6]    
+    f[5] = β * y[2] - μ * y[5]
+    f[6] = β * y[3] - μ * y[6]
     nothing
 end
 ps = (α₀=5e-4, α=.5, K=40.0, n=2, δ=(log(2)/120), β=(20*log(2)/120), μ=(log(2)/600))
@@ -77,10 +78,10 @@ i = indexof(os₁.P, states(oderepressilator))
 # @named rs₁ = ReactionSystem(rxs, t, specs, pars)
 # @named rs₂ = ReactionSystem(rxs, t, specs, pars)
 # @named rs₃ = ReactionSystem(rxs, t, specs, pars)
-# connections = [rs₁.R ~ rs₃.P, 
+# connections = [rs₁.R ~ rs₃.P,
 #                rs₂.R ~ rs₁.P,
 #                rs₃.R ~ rs₂.P]
 # @named csys = ODESystem(connections, t, [], [])
 # @named repressilator = ReactionSystem(t; systems=[csys,rs₁,rs₂,rs₃])
 # @named oderepressilator2 = convert(ODESystem, repressilator)
-# sys2 = structural_simplify(oderepressilator2)  # FAILS currently
+# sys2 = structural_simplify(oderepressilator2)  # FAILS currently