Merge branch 'master' into myb/alias

Author: YingboMa

src/systems/abstractsystem.jl

@@ -222,6 +222,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))
@@ -239,12 +250,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

@@ -192,39 +192,53 @@ function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
                                      sparse = false, simplify=false,
                                      kwargs...) where {iip}
 
-    idx = iip ? 2 : 1
-    f = generate_function(sys, dvs, ps; expression=Val{true}, kwargs...)[idx]
+    f_oop, f_iip = generate_function(sys, dvs, ps; expression=Val{true}, kwargs...)
+    fsym = gensym(:f)
+    _f = quote
+      $fsym(u,p,t) = $f_oop(u,p,t)
+      $fsym(du,u,p,t) = $f_iip(du,u,p,t)
+    end
+  
+    tgradsym = gensym(:tgrad)
     if tgrad
-        _tgrad = generate_tgrad(sys, dvs, ps;
+        tgrad_oop, tgrad_iip = generate_tgrad(sys, dvs, ps;
                                 simplify=simplify,
-                                expression=Val{true}, kwargs...)[idx]
+                                expression=Val{true}, kwargs...)
+        _tgrad = quote
+            $tgradsym(u,p,t) = $tgrad_oop(u,p,t)
+            $tgradsym(J,u,p,t) = $tgrad_iip(J,u,p,t)
+        end
     else
-        _tgrad = :nothing
+        _tgrad = :($tgradsym = nothing)
     end
 
+    jacsym = gensym(:jac)
     if jac
-        _jac = generate_jacobian(sys, dvs, ps;
+        jac_oop,jac_iip = generate_jacobian(sys, dvs, ps;
                                  sparse=sparse, simplify=simplify,
-                                 expression=Val{true}, kwargs...)[idx]
+                                 expression=Val{true}, kwargs...)
+        _jac = quote
+            $jacsym(u,p,t) = $jac_oop(u,p,t)
+            $jacsym(J,u,p,t) = $jac_iip(J,u,p,t)
+        end
     else
-        _jac = :nothing
+        _jac = :($jacsym = nothing)
     end
 
     M = calculate_massmatrix(sys)
 
     _M = (u0 === nothing || M == I) ? M : ArrayInterface.restructure(u0 .* u0',M)
 
     jp_expr = sparse ? :(similar($(get_jac(sys)[]),Float64)) : :nothing
-
     ex = quote
-        f = $f
-        tgrad = $_tgrad
-        jac = $_jac
+        $_f
+        $_tgrad
+        $_jac
         M = $_M
         ODEFunction{$iip}(
-                          f,
-                          jac = jac,
-                          tgrad = tgrad,
+                          $fsym,
+                          jac = $jacsym,
+                          tgrad = $tgradsym,
                           mass_matrix = M,
                           jac_prototype = $jp_expr,
                           syms = $(Symbol.(states(sys))),
@@ -252,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/diffeqs/odesystem.jl

@@ -224,3 +224,5 @@ function flatten(sys::ODESystem)
                         )
     end
 end
+
+ODESystem(eq::Equation, args...; kwargs...) = ODESystem([eq], args...; kwargs...)

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

@@ -35,21 +35,38 @@ generate_function(de, [x,y,z], [σ,ρ,β])
 jac_expr = generate_jacobian(de)
 jac = calculate_jacobian(de)
 jacfun = eval(jac_expr[2])
-# iip
-f = ODEFunction(de, [x,y,z], [σ,ρ,β])
-du = zeros(3)
-u  = collect(1:3)
-p  = collect(4:6)
-f(du, u, p, 0.1)
-@test du == [4, 0, -16]
-J = zeros(3, 3)
-jacfun(J, u, p, t)
-# oop
-f = ODEFunction(de, [x,y,z], [σ,ρ,β])
-du = @SArray zeros(3)
-u  = SVector(1:3...)
-p  = SVector(4:6...)
-@test f(u, p, 0.1) === @SArray [4, 0, -16]
+
+for f in [
+    ODEFunction(de, [x,y,z], [σ,ρ,β], tgrad = true, jac = true),
+    eval(ODEFunctionExpr(de, [x,y,z], [σ,ρ,β], tgrad = true, jac = true)),
+]
+    # iip
+    du = zeros(3)
+    u  = collect(1:3)
+    p  = collect(4:6)
+    f.f(du, u, p, 0.1)
+    @test du == [4, 0, -16]
+
+    # oop
+    du = @SArray zeros(3)
+    u  = SVector(1:3...)
+    p  = SVector(4:6...)
+    @test f.f(u, p, 0.1) === @SArray [4, 0, -16]
+
+    # iip vs oop
+    du = zeros(3)
+    g = similar(du)
+    J = zeros(3, 3)
+    u  = collect(1:3)
+    p  = collect(4:6)
+    f.f(du, u, p, 0.1)
+    @test du == f(u, p, 0.1)
+    f.tgrad(g, u, p, t)
+    @test g == f.tgrad(u, p, t)
+    f.jac(J, u, p, t)
+    @test J == f.jac(u, p, t)
+end
+
 
 eqs = [D(x) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y*t,
@@ -59,6 +76,8 @@ ModelingToolkit.calculate_tgrad(de)
 
 tgrad_oop, tgrad_iip = eval.(ModelingToolkit.generate_tgrad(de))
 
+u  = SVector(1:3...)
+p  = SVector(4:6...)
 @test tgrad_oop(u,p,t) == [0.0,-u[2],0.0]
 du = zeros(3)
 tgrad_iip(du,u,p,t)
@@ -254,3 +273,11 @@ sys = ODESystem(eqs, t)
 @test isequal(ModelingToolkit.get_iv(sys), t)
 @test isequal(states(sys), [x1, x2])
 @test isempty(parameters(sys))
+
+# one equation ODESystem test
+@parameters t r
+@variables x(t)
+D = Differential(t)
+eq = D(x) ~ r*x
+ode = ODESystem(eq)
+@test equations(ode) == [eq]

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())