diff --git a/src/reaction.jl b/src/reaction.jl index 23538214aa..701bbfebab 100644 --- a/src/reaction.jl +++ b/src/reaction.jl @@ -72,6 +72,14 @@ 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_vector = type[value(v) for v in vec] + return convert(Vector{type}, type_vector) +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 +93,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,7 +139,7 @@ Notes: - The three-argument form assumes all reactant and product stoichiometric coefficients are one. """ -struct Reaction{S, T} +struct Reaction{S,T} """The rate function (excluding mass action terms).""" rate::Any """Reaction substrates.""" @@ -160,83 +165,76 @@ 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{Q}, prods::Vector{R}, substoich::Vector{S}, prodstoich::Vector{T}; + netstoich = [], metadata = Pair{Symbol, Any}[], + only_use_rate = metadata_only_use_rate_check(metadata), kwargs...) where {Q,R,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. + reactant_type = promote_type(Q, R) + (reactant_type <: Num) && (reactant_type = Any) + 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{reactant_type, 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. diff --git a/test/reactionsystem_core/reaction.jl b/test/reactionsystem_core/reaction.jl index 4760641c4b..e70544d0fb 100644 --- a/test/reactionsystem_core/reaction.jl +++ b/test/reactionsystem_core/reaction.jl @@ -8,6 +8,89 @@ using ModelingToolkit: value, get_variables! # Sets the default `t` to use. t = default_t() +### Reaction Constructor Tests ### + +# Checks that `Reaction`s can be successfully 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) + + # Tries for different types of rates (should not matter). + for rate in (k, k*A, 2, 3.0, 4//3) + # Creates `Reaction`s. + rx1 = Reaction(rate, [X], []) + rx2 = Reaction(rate, [x], [Y], [1.5], [1]) + rx3 = Reaction(rate, [x, X], [], [n1 + n2, n2], []) + rx4 = Reaction(rate, [X, Y], [X, Y, Z], [2//3, 3], [1//3, 1, 2]) + rx5 = Reaction(rate, [X, Y], [X, Y, Z], [2, 3], [1, n1, n2]) + rx6 = Reaction(rate, [X], [x], [n1], [1]) + + # Check `Reaction` types. + @test rx1 isa Reaction{Any,Int64} + @test rx2 isa Reaction{Any,Float64} + @test rx3 isa Reaction{Any,Any} + @test rx4 isa Reaction{Any,Rational{Int64}} + @test rx5 isa Reaction{Any,Any} + @test rx6 isa Reaction{Any,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 +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 + + ### Test Basic Accessors ### # Tests the `get_variables` function. @@ -42,7 +125,6 @@ end # Tests basic accessor functions. # Tests that repeated metadata entries are not permitted. let - @variables t @parameters k @species X(t) X2(t) @@ -60,7 +142,6 @@ end # Tests accessors for system without metadata. let - @variables t @parameters k @species X(t) X2(t) @@ -77,7 +158,6 @@ end # Tests basic accessor functions. # Tests various metadata types. let - @variables t @parameters k @species X(t) X2(t) @@ -109,7 +189,6 @@ end # Tests the noise scaling metadata. let - @variables t @parameters k η @species X(t) X2(t)