Merge pull request #534 from isaacsas/fix_symbol_to_var

Fix symbol to variable mappings

Author: isaacsas

src/networkapi.jl

@@ -408,7 +408,7 @@ function setdefaults!(rn, newdefs)
     defs = eltype(newdefs) <: Pair{Symbol} ? symmap_to_varmap(rn, newdefs) : newdefs
     rndefs = MT.get_defaults(rn)
     for (var, val) in defs
-        rndefs[var] = value(val)
+        rndefs[value(var)] = value(val)
     end
     nothing
 end

src/reactionsystem.jl

@@ -690,6 +690,29 @@ function MT.observed(sys::ReactionSystem)
            init = obs)
 end
 
+function Base.propertynames(sys::ReactionSystem; private = false)
+    if private
+        return fieldnames(typeof(sys))
+    else
+        names = Symbol[]
+        for s in get_systems(sys)
+            push!(names, getname(s))
+        end
+        MT.has_states(sys) && for s in get_states(sys)
+            push!(names, getname(s))
+        end
+        MT.has_ps(sys) && for s in get_ps(sys)
+            push!(names, getname(s))
+        end
+        MT.has_observed(sys) && for s in get_observed(sys)
+            push!(names, getname(s.lhs))
+        end
+        (get_constraints(sys) !== nothing) &&
+            append!(names, propertynames(get_constraints(sys)))
+        return unique!(names)
+    end
+end
+
 function MT.getvar(sys::ReactionSystem, name::Symbol; namespace = false)
     if isdefined(sys, name)
         Base.depwarn("`sys.name` like `sys.$name` is deprecated. Use getters like `get_$name` instead.",
@@ -702,11 +725,6 @@ function MT.getvar(sys::ReactionSystem, name::Symbol; namespace = false)
     i === nothing ||
         return namespace ? rename(systems[i], renamespace(sys, name)) : systems[i]
 
-    constraints = get_constraints(sys)
-    if constraints !== nothing && nameof(constraints) == name
-        return namespace ? rename(constraints, renamespace(sys, name)) : constraints
-    end
-
     avs = get_var_to_name(sys)
     v = get(avs, name, nothing)
     v === nothing || return namespace ? renamespace(sys, v) : v
@@ -723,6 +741,7 @@ function MT.getvar(sys::ReactionSystem, name::Symbol; namespace = false)
     i = findfirst(x -> getname(x.lhs) == name, obs)
     i === nothing || return namespace ? renamespace(sys, obs[i]) : obs[i]
 
+    constraints = get_constraints(sys)
     if constraints !== nothing
         return getvar(rename(constraints, nameof(sys)), name; namespace = namespace)
     end

test/api.jl

@@ -586,6 +586,7 @@ setdefaults!(rn, [:S => 999.0, :I => 1.0, :R => 0.0, :α => 1e-4, :β => 0.01])
 op = ODEProblem(rn, [], tspan, [])
 sol2 = solve(op, Tsit5())
 @test norm(sol.u - sol2.u) ≈ 0
+@test all(p -> p[1] isa Symbolics.Symbolic, collect(ModelingToolkit.defaults(rn)))
 
 rn = @reaction_network begin
     α, S + I --> 2I

test/reactionsystem_components.jl

@@ -1,3 +1,4 @@
+#! format: off
 using ModelingToolkit, Catalyst, LinearAlgebra, OrdinaryDiffEq, Test
 using SciMLNLSolve
 
@@ -262,6 +263,13 @@ u₀ = [A => 1.0, B => 2.0, C => 0.0]
 oprob = ODEProblem(structural_simplify(osys), u₀, (0.0, 10.0), p)
 sol = solve(oprob, Tsit5())
 @test isapprox(0, norm(sol[ns.D] .- 2 * sol[A] - 3 * sol[B]), atol = (100 * eps()))
+psyms = [:r₊ => 1.0, :r₋ => 2.0, :ns₊β => 3.0]
+u₀syms = [:A => 1.0, :B => 2.0, :C => 0.0]
+p = symmap_to_varmap(osys, psyms)
+u₀ = symmap_to_varmap(osys, u₀syms)
+oprob = ODEProblem(structural_simplify(osys), u₀, (0.0, 10.0), p)
+sol = solve(oprob, Tsit5())
+@test isapprox(0, norm(sol[ns.D] .- 2 * sol[A] - 3 * sol[B]), atol = (100 * eps()))
 
 # test API functions for composed model
 @test issetequal(species(rs), [A, B, C])
@@ -354,3 +362,19 @@ rxs2 = [Reaction(p2, [ParentScope(A)], nothing)]
 rsc = compose(rs1, [rs2])
 orsc = convert(ODESystem, rsc)
 @test length(equations(orsc)) == 1
+
+# test constraint system symbols can be set via setdefaults!
+let
+    @parameters b
+    @variables t V(t) [isbcspecies = true]
+    rn = @reaction_network begin
+        k/$V, A + B --> C
+    end k
+    Dt = Differential(t)
+    @named csys = ODESystem([Dt(V) ~ -b * V], t)
+    @named fullrn = extend(csys, rn)
+    setdefaults!(fullrn, [:b => 2.0])
+    @unpack b = fullrn
+    @test haskey(ModelingToolkit.defaults(fullrn), b)
+    @test ModelingToolkit.defaults(fullrn)[b] == 2.0
+end