Merge pull request #546 from isaacsas/parametric_composition_bug

Parametric composition bug

Author: isaacsas

Project.toml

@@ -1,6 +1,6 @@
 name = "Catalyst"
 uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
-version = "12.2.0"
+version = "12.2.1"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"

src/reactionsystem.jl

@@ -187,14 +187,27 @@ function Base.show(io::IO, rx::Reaction)
     print_rxside(io, rx.products, rx.prodstoich)
 end
 
+function apply_if_nonempty(f, v)
+    isempty(v) && return v
+    s = similar(v)
+    map!(f, s, v)
+    s
+end
+
 function ModelingToolkit.namespace_equation(rx::Reaction, name)
-    subs = isempty(rx.substrates) ? rx.substrates :
-           [namespace_expr(sub, name) for sub in rx.substrates]
-    prods = isempty(rx.products) ? rx.products :
-            [namespace_expr(prod, name) for prod in rx.products]
-    Reaction(namespace_expr(rx.rate, name),
-             subs, prods, rx.substoich, rx.prodstoich,
-             [namespace_expr(n[1], name) => n[2] for n in rx.netstoich], rx.only_use_rate)
+    f = Base.Fix2(namespace_expr, name)
+    rate = f(rx.rate)
+    subs = apply_if_nonempty(f, rx.substrates)
+    prods = apply_if_nonempty(f, rx.products)
+    substoich = apply_if_nonempty(f, rx.substoich)
+    prodstoich = apply_if_nonempty(f, rx.prodstoich)
+    netstoich = if isempty(rx.netstoich)
+        rx.netstoich
+    else
+        ns = similar(rx.netstoich)
+        map!(n -> f(n[1]) => f(n[2]), ns, rx.netstoich)
+    end
+    Reaction(rate, subs, prods, substoich, prodstoich, netstoich, rx.only_use_rate)
 end
 
 netstoich_stoichtype(::Vector{Pair{S, T}}) where {S, T} = T

test/reactionsystem_components.jl

@@ -378,3 +378,25 @@ let
     @test haskey(ModelingToolkit.defaults(fullrn), b)
     @test ModelingToolkit.defaults(fullrn)[b] == 2.0
 end
+
+# https://github.com/SciML/Catalyst.jl/issues/545
+let
+    rn_AB = @reaction_network AB begin
+        k1, A --> n*B
+    end k1 n
+
+    rn_BC = @reaction_network BC begin
+        k2, B --> C
+    end k2
+
+    @variables t
+    @named rs = ReactionSystem(t; systems = [rn_AB, rn_BC])
+    sts = states(rs)
+    @test issetequal(sts, (@variables AB₊A(t) AB₊B(t) BC₊B(t) BC₊C(t)))
+    ps = parameters(rs)
+    @test issetequal(ps, (@parameters AB₊k1 AB₊n BC₊k2))
+    rxs = reactions(rs)
+    @parameters AB₊n
+    rxs2 = Reaction[(@reaction AB₊k1, AB₊A --> $(AB₊n)*AB₊B), (@reaction BC₊k2, BC₊B --> BC₊C)]
+    @test (length(rxs) == length(rxs2)) && issubset(rxs, rxs2)
+end