try to update scoping in constructor

Author: isaacsas

src/reaction.jl

@@ -351,20 +351,16 @@ MT.is_alg_equation(rx::Reaction) = false
 # MTK functions for extracting variables within equation type object
 MT.eqtype_supports_collect_vars(rx::Reaction) = true
 function MT.collect_vars!(unknowns, parameters, rx::Reaction, iv; depth = 0,
-        op = MT.Differential)
+        op = MT.Differential)    
+
     MT.collect_vars!(unknowns, parameters, rx.rate, iv; depth, op)
-    for sub in rx.substrates
-        MT.collect_vars!(unknowns, parameters, sub, iv; depth, op)
-    end
-    for prod in rx.products
-        MT.collect_vars!(unknowns, parameters, prod, iv; depth, op)
-    end
-    for substoich in rx.substoich
-        MT.collect_vars!(unknowns, parameters, substoich, iv; depth, op)
-    end
-    for prodstoich in rx.prodstoich
-        MT.collect_vars!(unknowns, parameters, prodstoich, iv; depth, op)
+
+    for items in (rx.substrates, rx.products, rx.substoich, rx.prodstoich)
+        for item in items
+            MT.collect_vars!(unknowns, parameters, item, iv; depth, op)
+        end
     end
+
     if hasnoisescaling(rx)
         ns = getnoisescaling(rx)
         MT.collect_vars!(unknowns, parameters, ns, iv; depth, op)

src/reactionsystem.jl

@@ -512,24 +512,8 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in;
     eqs = Equation[eq for eq in rxs_and_eqs if eq isa Equation]
 
     # Loops through all reactions, adding encountered quantities to the unknown and parameter vectors.
-    # Starts by looping through substrates + products only (so these are added to the vector first).
-    # Next, the other components of reactions (e.g. rates and stoichiometries) are added.
     for rx in rxs
-        for reactants in (rx.substrates, rx.products), spec in reactants
-            MT.isparameter(spec) ? push!(ps, spec) : push!(us, spec)
-        end
-    end
-    for rx in rxs
-        # Adds all quantities encountered in the reaction's rate.
-        findvars!(ps, us, rx.rate, ivs, vars)
-
-        # Extracts all quantities encountered within stoichiometries.
-        for stoichiometry in (rx.substoich, rx.prodstoich), sym in stoichiometry
-            (sym isa Symbolic) && findvars!(ps, us, sym, ivs, vars)
-        end
-
-        # Extract all quantities encountered in relevant `Reaction` metadata.
-        hasnoisescaling(rx) && findvars!(ps, us, getnoisescaling(rx), ivs, vars)
+        MT.collect_vars!(us, ps, rx, iv)
     end
 
     # Extracts any species, variables, and parameters that occur in (non-reaction) equations.
@@ -543,6 +527,11 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in;
         fulleqs = rxs
     end
 
+    # get variables in subsystems with scope at this level
+    for ssys in get(kwargs, :systems, [])
+        MT.collect_scoped_vars!(us, ps, ssys, iv)
+    end
+
     # Loops through all events, adding encountered quantities to the unknown and parameter vectors.
     find_event_vars!(ps, us, continuous_events, ivs, vars)
     find_event_vars!(ps, us, discrete_events, ivs, vars)

test/compositional_modelling/component_based_model_creation.jl

@@ -501,3 +501,46 @@ let
     defs[rn2.X] == 30.0
     defs[rn2.Z] == 40.0
 end
+
+# test scoping in compose
+# code adapted from ModelingToolkit.jl tests
+let
+    t = default_t()
+    D = default_time_deriv()
+    @species x1(t) x2(t) 
+    @variables x3(t) x4(t) x5(t)
+    x2 = ParentScope(x2)
+    x3 = ParentScope(ParentScope(x3))
+    x4 = DelayParentScope(x4, 2)
+    x5 = GlobalScope(x5)
+    @parameters p1 p2 p3 p4 p5
+    p2 = ParentScope(p2)
+    p3 = ParentScope(ParentScope(p3))
+    p4 = DelayParentScope(p4, 2)
+    p5 = GlobalScope(p5)
+    rxs = [Reaction(p1, nothing, [x1]), Reaction(p2, [x2], nothing), 
+           D(x3) ~ p3, D(x4) ~ p4, D(x5) ~ p5]
+    @named sys1 = ReactionSystem(rxs, t)
+    @test isequal(x1, only(unknowns(sys1)))
+    @test isequal(x1, only(species(sys1)))
+    @test isequal(p1, only(parameters(sys1)))
+    @named sys2 = ReactionSystem([], t; systems = [sys1])
+    @test length(unknowns(sys2)) == 2
+    @test any(isequal(x2), unknowns(sys2))
+    @test any(isequal(x2), species(sys2))
+    @test length(parameters(sys2)) == 2
+    @test any(isequal(p2), parameters(sys2))
+    @named sys3 = ReactionSystem(Equation[], t)
+    sys3 = sys3 ∘ sys2
+    @test length(unknowns(sys3)) == 4
+    @test any(isequal(x3), unknowns(sys3))
+    @test any(isequal(x4), unknowns(sys3))
+    @test length(species(sys3)) == 2
+    @test length(parameters(sys3)) == 4
+    @test any(isequal(p3), parameters(sys3))
+    @test any(isequal(p4), parameters(sys3))
+    sys4 = complete(sys3)
+    @test length(unknowns(sys3)) == 4
+    @test length(parameters(sys4)) == 5
+    @test any(isequal(p5), parameters(sys4))
+end