Merge pull request #856 from SciML/update_possible_Reaction_inputs

Update `Reaction` constructor

Author: TorkelE

src/reaction.jl

@@ -72,6 +72,13 @@ 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
@@ -85,9 +92,6 @@ 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 ###
 
 """
@@ -134,19 +138,19 @@ Notes:
 - The three-argument form assumes all reactant and product stoichiometric coefficients
   are one.
 """
-struct Reaction{S, T}
+struct Reaction{T}
     """The rate function (excluding mass action terms)."""
     rate::Any
     """Reaction substrates."""
-    substrates::Vector
+    substrates::Vector{Any}
     """Reaction products."""
-    products::Vector
+    products::Vector{Any}
     """The stoichiometric coefficients of the reactants."""
-    substoich::Vector{T}
+    substoich::Vector
     """The stoichiometric coefficients of the products."""
-    prodstoich::Vector{T}
+    prodstoich::Vector
     """The net stoichiometric coefficients of all species changed by the reaction."""
-    netstoich::Vector{Pair{S, T}}
+    netstoich::Vector{Pair{Any, 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.
@@ -160,83 +164,74 @@ struct Reaction{S, T}
 end
 
 # Five-argument constructor accepting rate, substrates, and products, and their stoichiometries.
-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
+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."))
     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)
 
-    # try to get a common type for stoichiometry, using Any if have Syms
+    # Ensures everything have uniform and correct types.
+    subs = promote_reaction_vector(subs, BasicSymbolic{Real})
+    prods =  promote_reaction_vector(prods, BasicSymbolic{Real})
     stoich_type = promote_type(S, T)
-    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
+    (stoich_type <: Num) && (stoich_type = Any)
+    substoich = promote_reaction_vector(substoich, stoich_type)
+    prodstoich = promote_reaction_vector(prodstoich, stoich_type)
 
+    # 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
-
-    ns = if netstoich === nothing
-        get_netstoich(subs, prods, substoich′, prodstoich′)
-    else
-        (netstoich_stoichtype(netstoich) != stoich_type) ?
-        convert.(stoich_type, netstoich) : netstoich
+    
+    # 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]
     end
 
-    # Check that all metadata entries are unique. (cannot use `in` since some entries may be symbolics).
+    # Handles metadata (check that all entries are unique, remove potential `only_use_rate`
+    # entries, and converts to the required type.
     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
-
-    # Ensures metadata have the correct type.
     metadata = convert(Vector{Pair{Symbol, Any}}, metadata)
 
-    Reaction(value(rate), subs, prods, substoich′, prodstoich′, ns, only_use_rate, metadata)
+    Reaction{stoich_type}(value(rate), subs, prods, substoich, prodstoich, netstoich, only_use_rate, metadata)
 end
 
-# 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...)
+# 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...)
+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...)
 end
 
 # Union type for `Reaction`s and `Equation`s.

test/reactionsystem_core/reaction.jl

@@ -36,13 +36,94 @@ 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)
 
@@ -60,7 +141,6 @@ end
 
 # Tests accessors for system without metadata.
 let
-    @variables t
     @parameters k
     @species X(t) X2(t)
 
@@ -77,7 +157,6 @@ end
 # Tests basic accessor functions.
 # Tests various metadata types.
 let
-    @variables t
     @parameters k
     @species X(t) X2(t)
 
@@ -109,7 +188,6 @@ end
 
 # Tests the noise scaling metadata.
 let
-    @variables t
     @parameters k η  
     @species X(t) X2(t)
 

test/reactionsystem_core/reactionsystem.jl

@@ -140,7 +140,6 @@ 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)