Revert "Update Reaction constructor"

SciML · Jun 5, 2024 · ea99c8b · ea99c8b
1 parent ae6cf0a
commit ea99c8b
Show file tree

Hide file tree

Showing 3 changed files with 69 additions and 141 deletions.
diff --git a/src/reaction.jl b/src/reaction.jl
@@ -72,13 +72,6 @@ function metadata_only_use_rate_check(metadata)
     return Bool(metadata[only_use_rate_idx][2])
 end
 
-# Used to promote a vector to the appropriate type. Takes the `vec == Any[]` case into account by
-# returning an empty vector of the appropriate type.
-function promote_reaction_vector(vec, type)
-    isempty(vec) && (return type[])
-    type[value(v) for v in vec]
-end
-
 # calculates the net stoichiometry of a reaction as a vector of pairs (sub,substoich)
 function get_netstoich(subs, prods, sstoich, pstoich)
     # stoichiometry as a Dictionary
@@ -92,6 +85,9 @@ function get_netstoich(subs, prods, sstoich, pstoich)
     [el for el in nsdict if !_iszero(el[2])]
 end
 
+# Get the net stoichiometries' type.
+netstoich_stoichtype(::Vector{Pair{S, T}}) where {S, T} = T
+
 ### Reaction Structure ###
 
 """
@@ -138,19 +134,19 @@ Notes:
 - The three-argument form assumes all reactant and product stoichiometric coefficients
   are one.
 """
-struct Reaction{T}
+struct Reaction{S, T}
     """The rate function (excluding mass action terms)."""
     rate::Any
     """Reaction substrates."""
-    substrates::Vector{Any}
+    substrates::Vector
     """Reaction products."""
-    products::Vector{Any}
+    products::Vector
     """The stoichiometric coefficients of the reactants."""
-    substoich::Vector
+    substoich::Vector{T}
     """The stoichiometric coefficients of the products."""
-    prodstoich::Vector
+    prodstoich::Vector{T}
     """The net stoichiometric coefficients of all species changed by the reaction."""
-    netstoich::Vector{Pair{Any, T}}
+    netstoich::Vector{Pair{S, T}}
     """
     `false` (default) if `rate` should be multiplied by mass action terms to give the rate law.
     `true` if `rate` represents the full reaction rate law.
@@ -164,74 +160,83 @@ struct Reaction{T}
 end
 
 # Five-argument constructor accepting rate, substrates, and products, and their stoichiometries.
-function Reaction(rate, subs::Vector, prods::Vector, substoich::Vector{S}, prodstoich::Vector{T};
-                   netstoich = [], metadata = Pair{Symbol, Any}[], 
-                   only_use_rate = metadata_only_use_rate_check(metadata), kwargs...) where {S,T}
-    # Error checks.
-    isempty(subs) && isempty(prods) &&
-        throw(ArgumentError("A reaction requires either a non-empty substrate or product vector."))
-    length(subs) != length(substoich) &&
-        throw(ArgumentError("The substrate vector ($(subs)) and the substrate stoichiometry vector ($(substoich)) must have equal length."))
-    length(prods) != length(prodstoich) &&
-        throw(ArgumentError("The product vector ($(prods)) and the product stoichiometry vector ($(prodstoich)) must have equal length."))
+function Reaction(rate, subs, prods, substoich, prodstoich;
+                  netstoich = nothing, metadata = Pair{Symbol, Any}[], 
+                  only_use_rate = metadata_only_use_rate_check(metadata), kwargs...)
+    (isnothing(prods) && isnothing(subs)) &&
+        throw(ArgumentError("A reaction requires a non-nothing substrate or product vector."))
+    (isnothing(prodstoich) && isnothing(substoich)) &&
+        throw(ArgumentError("Both substrate and product stochiometry inputs cannot be nothing."))
+
+    if isnothing(subs)
+        prodtype = typeof(value(first(prods)))
+        subs = Vector{prodtype}()
+        !isnothing(substoich) &&
+            throw(ArgumentError("If substrates are nothing, substrate stoichiometries have to be so too."))
+        substoich = typeof(prodstoich)()
+    else
+        subs = value.(subs)
+    end
     allunique(subs) ||
         throw(ArgumentError("Substrates can not be repeated in the list provided to `Reaction`, please modify the stoichiometry for any repeated substrates instead."))
+    S = eltype(substoich)
+
+    if isnothing(prods)
+        prods = Vector{eltype(subs)}()
+        !isnothing(prodstoich) &&
+            throw(ArgumentError("If products are nothing, product stoichiometries have to be so too."))
+        prodstoich = typeof(substoich)()
+    else
+        prods = value.(prods)
+    end
     allunique(prods) ||
         throw(ArgumentError("Products can not be repeated in the list provided to `Reaction`, please modify the stoichiometry for any repeated products instead."))
+    T = eltype(prodstoich)
 
-    # Ensures everything have uniform and correct types.
-    subs = promote_reaction_vector(subs, BasicSymbolic{Real})
-    prods =  promote_reaction_vector(prods, BasicSymbolic{Real})
+    # try to get a common type for stoichiometry, using Any if have Syms
     stoich_type = promote_type(S, T)
-    (stoich_type <: Num) && (stoich_type = Any)
-    substoich = promote_reaction_vector(substoich, stoich_type)
-    prodstoich = promote_reaction_vector(prodstoich, stoich_type)
+    if stoich_type <: Num
+        stoich_type = Any
+        substoich′ = Any[value(s) for s in substoich]
+        prodstoich′ = Any[value(p) for p in prodstoich]
+    else
+        substoich′ = (S == stoich_type) ? substoich : convert.(stoich_type, substoich)
+        prodstoich′ = (T == stoich_type) ? prodstoich : convert.(stoich_type, prodstoich)
+    end
 
-    # Checks that all reactants are valid.
     if !(all(isvalidreactant, subs) && all(isvalidreactant, prods))
         badsts = union(filter(!isvalidreactant, subs), filter(!isvalidreactant, prods))
-        throw(ArgumentError("To be a valid substrate or product, non-constant species must be declared via @species, while constant species must be parameters with the isconstantspecies metadata. The following reactants do not follow this convention:\n $badsts"))
+        throw(ArgumentError("""To be a valid substrate or product, non-constant species must be declared via @species, while constant species must be parameters with the isconstantspecies metadata. The following reactants do not follow this convention:\n $badsts"""))
     end
-
-    # Computes the net stoichiometries.
-    if isempty(netstoich)
-        netstoich = get_netstoich(subs, prods, substoich, prodstoich)
-    elseif typeof(netstoich) != Vector{Pair{BasicSymbolic{Real}, stoich_type}}
-        netstoich = Pair{BasicSymbolic{Real}, stoich_type}[
-                        value(ns[1]) => convert(stoich_type, ns[2]) for ns in netstoich]
+
+    ns = if netstoich === nothing
+        get_netstoich(subs, prods, substoich′, prodstoich′)
+    else
+        (netstoich_stoichtype(netstoich) != stoich_type) ?
+        convert.(stoich_type, netstoich) : netstoich
     end
 
-    # Handles metadata (check that all entries are unique, remove potential `only_use_rate`
-    # entries, and converts to the required type.
+    # Check that all metadata entries are unique. (cannot use `in` since some entries may be symbolics).
     if !allunique(entry[1] for entry in metadata)
         error("Repeated entries for the same metadata encountered in the following metadata set: $([entry[1] for entry in metadata]).")
     end
+
+    # Deletes potential `:only_use_rate => ` entries from the metadata.
     if any(:only_use_rate == entry[1] for entry in metadata) 
         deleteat!(metadata, findfirst(:only_use_rate == entry[1] for entry in metadata))
     end
-    metadata = convert(Vector{Pair{Symbol, Any}}, metadata)
 
-    Reaction{stoich_type}(value(rate), subs, prods, substoich, prodstoich, netstoich, only_use_rate, metadata)
-end
+    # Ensures metadata have the correct type.
+    metadata = convert(Vector{Pair{Symbol, Any}}, metadata)
 
-# Three-argument constructor. Handles the case where no stoichiometries is given
-# (by assuming that all stoichiometries are `1`).
-function Reaction(rate, subs::Vector, prods::Vector; kwargs...)
-    Reaction(rate, subs, prods, ones(Int, length(subs)), ones(Int, length(prods)); kwargs...)
+    Reaction(value(rate), subs, prods, substoich′, prodstoich′, ns, only_use_rate, metadata)
 end
 
-# Handles cases where `nothing` is given (instead of an empty vector). 
-function Reaction(rate, subs::Vector, prods::Nothing, substoich::Vector, prodstoich::Nothing; kwargs...)
-    Reaction(rate, subs, Int64[], substoich, Int64[]; kwargs...)
-end
-function Reaction(rate, subs::Nothing, prods::Vector, substoich::Nothing, prodstoich::Vector; kwargs...)
-    Reaction(rate, Int64[], prods, Int64[], prodstoich; kwargs...)
-end
-function Reaction(rate, subs::Vector, prods::Nothing; kwargs...)
-    Reaction(rate, subs, Int64[]; kwargs...)
-end
-function Reaction(rate, subs::Nothing, prods::Vector; kwargs...)
-    Reaction(rate, Int64[], prods; kwargs...)
+# Three argument constructor assumes stoichiometric coefs are one and integers.
+function Reaction(rate, subs, prods; kwargs...)
+    sstoich = isnothing(subs) ? nothing : ones(Int, length(subs))
+    pstoich = isnothing(prods) ? nothing : ones(Int, length(prods))
+    Reaction(rate, subs, prods, sstoich, pstoich; kwargs...)
 end
 
 # Union type for `Reaction`s and `Equation`s.

diff --git a/test/reactionsystem_core/reaction.jl b/test/reactionsystem_core/reaction.jl
@@ -36,94 +36,13 @@ let
     @test issetequal(get_symbolics(rx4), [X, t, Y, η1, η2])
 end
 
-### Reaction Constructor Tests ###
-
-# Sets the default `t` to use.
-t = default_t()
-
-# Checks that `Reaction`s can be sucesfully created using various complicated inputs.
-# Checks that the `Reaction`s have the correct type, and the correct net stoichiometries are generated.
-let
-    # Declare symbolic variables.
-    @parameters k n1 n2::Int32 x [isconstantspecies=true]
-    @species X(t) Y(t) Z(t)
-    @variables A(t)
-
-    # Creates `Reaction`s.
-    rx1 = Reaction(k*A, [X], [])
-    rx2 = Reaction(k*A, [x], [Y], [1.5], [1])
-    rx3 = Reaction(k*A, [x, X], [], [n1 + n2, n2], [])
-    rx4 = Reaction(k*A, [X, Y], [X, Y, Z], [2//3, 3], [1//3, 1, 2])
-    rx5 = Reaction(k*A, [X, Y], [X, Y, Z], [2, 3], [1, n1, n2])
-    rx6 = Reaction(k*A, [X], [x], [n1], [1])
-
-    # Check `Reaction` types.
-    @test rx1 isa Reaction{Int64}
-    @test rx2 isa Reaction{Float64}
-    @test rx3 isa Reaction{Any}
-    @test rx4 isa Reaction{Rational{Int64}}
-    @test rx5 isa Reaction{Any}
-    @test rx6 isa Reaction{Any}
-
-    # Check `Reaction` net stoichiometries.
-    issetequal(rx1.netstoich, [X => -1])
-    issetequal(rx2.netstoich, [x => -1.5, Y => 1.0])
-    issetequal(rx3.netstoich, [x => -n1 - n2, X => -n2])
-    issetequal(rx4.netstoich, [X => -1//3, Y => -2//1, Z => 2//1])
-    issetequal(rx5.netstoich, [X => -1, Y => n1 - 3, Z => n2])
-    issetequal(rx6.netstoich, [X => -n1, x => 1])
-end
-
-# Tests that various `Reaction` constructors gives identical inputs.
-let
-    # Declare symbolic variables.
-    @parameters k n1 n2::Int32
-    @species X(t) Y(t) Z(t)
-    @variables A(t)
-
-    # Tests that the three-argument constructor generates correct result.
-    @test Reaction(k*A, [X], [Y, Z]) == Reaction(k*A, [X], [Y, Z], [1], [1, 1])
-
-    # Tests that `[]` and `nothing` can be used interchangeably.
-    @test Reaction(k*A, [X, Z], nothing) == Reaction(k*A, [X, Z], [])
-    @test Reaction(k*A, nothing, [Y, Z]) == Reaction(k*A, [], [Y, Z])
-    @test Reaction(k*A, [X, Z], nothing, [n1 + n2, 2], nothing) == Reaction(k*A, [X, Z], [], [n1 + n2, 2], [])
-    @test Reaction(k*A, nothing, [Y, Z], nothing, [n1 + n2, 2]) == Reaction(k*A, [], [Y, Z], [], [n1 + n2, 2])
-end
-
-# Tests that various incorrect inputs yields errors.
-let
-    # Declare symbolic variables.
-    @parameters k n1 n2::Int32
-    @species X(t) Y(t) Z(t)
-    @variables A(t)
-
-    # Neither substrates nor products.
-    @test_throws ArgumentError Reaction(k*A, [], [])
-
-    # Substrate vector not of equal length to substrate stoichiometry vector.
-    @test_throws ArgumentError Reaction(k*A, [X, X, Z], [], [1, 2], [])
-
-    # Product vector not of equal length to product stoichiometry vector.
-    @test_throws ArgumentError Reaction(k*A, [], [X, X, Z], [], [1, 2])
-
-    # Repeated substrates.
-    @test_throws ArgumentError Reaction(k*A, [X, X, Z], [])
-
-    # Repeated products.
-    @test_throws ArgumentError Reaction(k*A, [], [Y, Z, Z])
-
-    # Non-valid reactants (parameter or variable).
-    @test_throws ArgumentError Reaction(k*A, [], [A])
-    @test_throws ArgumentError Reaction(k*A, [], [k])
-end
-
 ### Tests Metadata ###
 
 # Tests creation.
 # Tests basic accessor functions.
 # Tests that repeated metadata entries are not permitted.
 let
+    @variables t
     @parameters k
     @species X(t) X2(t)
 
@@ -141,6 +60,7 @@ end
 
 # Tests accessors for system without metadata.
 let
+    @variables t
     @parameters k
     @species X(t) X2(t)
 
@@ -157,6 +77,7 @@ end
 # Tests basic accessor functions.
 # Tests various metadata types.
 let
+    @variables t
     @parameters k
     @species X(t) X2(t)
 
@@ -188,6 +109,7 @@ end
 
 # Tests the noise scaling metadata.
 let
+    @variables t
     @parameters k η  
     @species X(t) X2(t)
 

diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl
@@ -140,6 +140,7 @@ end
 ### Check ODE, SDE, and Jump Functions ###
 
 # Test by evaluating drift and diffusion terms.
+
 let
     u = rnd_u0(rs, rng)
     p = rnd_ps(rs, rng)