Add get_variables dispatch for Reaction

src/reaction.jl

@@ -343,6 +343,45 @@
 MT.is_diff_equation(rx::Reaction) = false
 MT.is_alg_equation(rx::Reaction) = false
 
+"""
+    get_symbolics(set, rx::Reaction)
+
+Returns all symbolic variables that are part of a reaction. This includes all variables 
+encountered in:
+    - Rates.
+    - Among substrates and products.
+    - Among stoichiometries.
+    - Among potential noise scaling metadata.
+"""
+function get_symbolics(rx::Reaction)
+    return ModelingToolkit.get_variables!([], rx)
+end
+
+"""
+    get_variables!(set, rx::Reaction)
+
+Adds all symbolic variables that are part of a reaction to set. This includes all variables 
+encountered in:
+    - Rates.
+    - Among substrates and products.
+    - Among stoichiometries.
+    - Among potential noise scaling metadata.
+"""
+function ModelingToolkit.get_variables!(set, rx::Reaction)
+    if isdefined(Main, :Infiltrator)
+        Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__)
+      end
+    get_variables!(set, rx.rate)
+    foreach(sub -> push!(set, sub), rx.substrates)
+    foreach(prod -> push!(set, prod), rx.products)
+    for stoichs in (rx.substoich, rx.prodstoich), stoich in stoichs
+        (stoich isa BasicSymbolic) && get_variables!(set, stoich)
+    end
+    if has_noise_scaling(rx)
+        get_variables!(set, get_noise_scaling(rx))
+    end
+    return (set isa AbstractVector) ? unique!(set) : set
+end
 
 ### Dependency-related Functions ###
 

test/reactionsystem_core/reaction.jl

@@ -2,6 +2,39 @@
 
 # Fetch packages.
 using Catalyst, Test
+using Catalyst: get_symbolics
+using ModelingToolkit: value, get_variables!
+
+# Sets the default `t` to use.
+t = default_t()
+
+### Test Basic Accessors ###
+
+# Tests the `get_variables` function.
+let
+    # Declare symbolic variables.
+    @parameters k1 k2 n1 n2 η1 η2 p
+    @species X(t) Y(t) Z(t)
+    @variables A(t)
+
+    # Create `Reaction`s.
+    rx1 = Reaction(k1, [], [X])
+    rx2 = Reaction(k1 + k2, [X], [Y], [1], [n1]; metadata = [:noise_scaling => η1])
+    rx3 = Reaction(k1 + k2 + A, [X], [X, Y, Z], [1], [n1 + n2, 2, 1])
+    rx4 = Reaction(X + t, [], [Y]; metadata = [:noise_scaling => η1 + η2])
+
+    # Test `get_variables!`.
+    @test issetequal(get_variables!([value(p)], rx1), [k1, X, p])
+    @test issetequal(get_variables!([value(p)], rx2), [k1, k2, X, Y, n1, η1, p])
+    @test issetequal(get_variables!([value(p)], rx3), [k1, k2, A, X, Y, Z, n1, n2, p])
+    @test issetequal(get_variables!([value(p)], rx4), [X, t, Y, η1, η2, p])
+
+    # Test `get_symbolics`.
+    @test issetequal(get_symbolics(rx1), [k1, X])
+    @test issetequal(get_symbolics(rx2), [k1, k2, X, Y, n1, η1])
+    @test issetequal(get_symbolics(rx3), [k1, k2, A, X, Y, Z, n1, n2])
+    @test issetequal(get_symbolics(rx4), [X, t, Y, η1, η2])
+end
 
 ### Tests Metadata ###