allow setdefaults on constraint sys variables

Author: isaacsas

src/reactionsystem.jl

@@ -690,6 +690,28 @@ 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 +724,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 +740,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/reactionsystem_components.jl

@@ -262,6 +262,14 @@ 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