Merge pull request #777 from pepijndevos/sympar

Add support for symbolic parameters

Author: ChrisRackauckas

src/systems/abstractsystem.jl

@@ -221,6 +221,17 @@ function Base.getproperty(sys::AbstractSystem, name::Symbol)
     throw(error("Variable $name does not exist"))
 end
 
+function Base.setproperty!(sys::AbstractSystem, prop::Symbol, val)
+    if (pa = Sym{Parameter{Real}}(prop); pa in parameters(sys))
+        sys.default_p[pa] = value(val)
+    # comparing a Sym returns a symbolic expression
+    elseif (st = Sym{Real}(prop); any(s->s.name==st.name, states(sys)))
+        sys.default_u0[st] = value(val)
+    else
+        setfield!(sys, prop, val)
+    end
+end
+
 function renamespace(namespace, x)
     if x isa Num
         renamespace(namespace, value(x))
@@ -238,12 +249,12 @@ namespace_parameters(sys::AbstractSystem) = parameters(sys, parameters(sys))
 
 function namespace_default_u0(sys)
     d_u0 = default_u0(sys)
-    Dict(states(sys, k) => d_u0[k] for k in keys(d_u0))
+    Dict(states(sys, k) => namespace_expr(d_u0[k], nameof(sys), independent_variable(sys)) for k in keys(d_u0))
 end
 
 function namespace_default_p(sys)
     d_p = default_p(sys)
-    Dict(parameters(sys, k) => d_p[k] for k in keys(d_p))
+    Dict(parameters(sys, k) => namespace_expr(d_p[k], nameof(sys), independent_variable(sys)) for k in keys(d_p))
 end
 
 function namespace_equations(sys::AbstractSystem)

src/systems/diffeqs/abstractodesystem.jl

@@ -266,9 +266,12 @@ function process_DEProblem(constructor, sys::AbstractODESystem,u0map,parammap;
         u0 = nothing
     end
 
+    defp = default_p(sys)
     if !(parammap isa DiffEqBase.NullParameters)
         parammap′ = lower_mapnames(parammap)
-        p = varmap_to_vars(parammap′,ps; defaults=default_p(sys))
+        p = varmap_to_vars(parammap′,ps; defaults=defp)
+    elseif !isempty(defp)
+        p = varmap_to_vars(Dict(),ps; defaults=defp)
     else
         p = ps
     end

src/systems/nonlinear/nonlinearsystem.jl

@@ -211,10 +211,13 @@ function process_NonlinearProblem(constructor, sys::NonlinearSystem,u0map,paramm
     ps = parameters(sys)
     u0map′ = lower_mapnames(u0map)
     u0 = varmap_to_vars(u0map′,dvs; defaults=default_u0(sys))
+    defp = default_p(sys)
 
     if !(parammap isa DiffEqBase.NullParameters)
         parammap′ = lower_mapnames(parammap)
-        p = varmap_to_vars(parammap′,ps; defaults=default_p(sys))
+        p = varmap_to_vars(parammap′,ps; defaults=defp)
+    elseif !isempty(defp)
+        p = varmap_to_vars(Dict(),ps; defaults=defp)
     else
         p = ps
     end

src/variables.jl

@@ -213,6 +213,10 @@ applicable.
 function varmap_to_vars(varmap::Dict, varlist; defaults=Dict())
     varmap = merge(defaults, varmap) # prefers the `varmap`
     varmap = Dict(value(k)=>value(varmap[k]) for k in keys(varmap))
+    # resolve symbolic parameter expressions
+    for (p, v) in pairs(varmap)
+        varmap[p] = fixpoint_sub(v, varmap)
+    end
     T′ = eltype(values(varmap))
     T = Base.isconcretetype(T′) ? T′ : Base.promote_typeof(values(varmap)...)
     out = Vector{T}(undef, length(varlist))

test/runtests.jl

@@ -1,5 +1,6 @@
 using SafeTestsets, Test
 
+@safetestset "Symbolic parameters test" begin include("symbolic_parameters.jl") end
 @safetestset "Parsing Test" begin include("variable_parsing.jl") end
 @safetestset "Differentiation Test" begin include("derivatives.jl") end
 @safetestset "Simplify Test" begin include("simplify.jl") end

test/symbolic_parameters.jl

@@ -0,0 +1,48 @@
+using ModelingToolkit
+using NonlinearSolve
+using Test
+
+@variables x y z u
+@parameters σ ρ β
+
+eqs = [0 ~ σ*(y-x),
+       0 ~ x*(ρ-z)-y,
+       0 ~ x*y - β*z]
+
+par = [
+    σ => 1,
+    ρ => 0.1+σ,
+    β => ρ*1.1
+]
+u0 = Pair{Num, Any}[
+    x => u,
+    y => u,
+    z => u-0.1,
+]
+ns = NonlinearSystem(eqs, [x,y,z],[σ,ρ,β], name=:ns, default_p=par, default_u0=u0)
+ns.y = u*1.1
+ModelingToolkit.default_p(ns)
+resolved = ModelingToolkit.varmap_to_vars(Dict(), parameters(ns), defaults=ModelingToolkit.default_p(ns))
+@test resolved == [1, 0.1+1, (0.1+1)*1.1]
+
+prob = NonlinearProblem(ns, [u=>1.0], Pair[])
+@test prob.u0 == [1.0, 1.1, 0.9]
+@show sol = solve(prob,NewtonRaphson())
+
+@variables a
+@parameters b
+top = NonlinearSystem([0 ~ -a + ns.x+b], [a], [b], systems=[ns], name=:top)
+top.b = ns.σ*0.5
+top.ns.x = u*0.5
+
+res = ModelingToolkit.varmap_to_vars(Dict(), parameters(top), defaults=ModelingToolkit.default_p(top))
+@test res == [0.5, 1, 0.1+1, (0.1+1)*1.1]
+
+prob = NonlinearProblem(top, [states(ns, u)=>1.0, a=>1.0], Pair[])
+@test prob.u0 == [1.0, 0.5, 1.1, 0.9]
+@show sol = solve(prob,NewtonRaphson())
+
+# test NullParameters+defaults
+prob = NonlinearProblem(top, [states(ns, u)=>1.0, a=>1.0])
+@test prob.u0 == [1.0, 0.5, 1.1, 0.9]
+@show sol = solve(prob,NewtonRaphson())