Support for Jump & Reaction systems.

Author: lamorton

src/ModelingToolkit.jl

@@ -116,7 +116,6 @@ include("systems/diffeqs/sdesystem.jl")
 include("systems/diffeqs/abstractodesystem.jl")
 include("systems/diffeqs/first_order_transform.jl")
 include("systems/diffeqs/modelingtoolkitize.jl")
-include("systems/diffeqs/validation.jl")
 include("systems/diffeqs/basic_transformations.jl")
 
 include("systems/jumps/jumpsystem.jl")
@@ -130,10 +129,10 @@ include("systems/control/controlsystem.jl")
 include("systems/pde/pdesystem.jl")
 
 include("systems/reaction/reactionsystem.jl")
-include("systems/dependency_graphs.jl")
 
 include("systems/discrete_system/discrete_system.jl")
-
+include("systems/validation.jl")
+include("systems/dependency_graphs.jl")
 include("systems/systemstructure.jl")
 using .SystemStructures
 

src/systems/jumps/jumpsystem.jl

@@ -57,6 +57,7 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractSystem
         if checks
             check_variables(states, iv)
             check_parameters(ps, iv)
+            check_units(ap,iv)
         end
         new{U}(ap, iv, states, ps, var_to_name, observed, name, systems, defaults, connection_type)
     end

src/systems/reaction/reactionsystem.jl

@@ -162,6 +162,7 @@ struct ReactionSystem <: AbstractSystem
             ps′ = value.(ps)
             check_variables(states′, iv′)
             check_parameters(ps′, iv′)
+            check_units(eqs)
         end
         new(collect(eqs), iv′, states′, ps′, observed, name, systems, defaults, connection_type)
     end

src/systems/diffeqs/validation.jl renamed to src/systems/validation.jl

@@ -1,7 +1,7 @@
 Base.:*(x::Union{Num,Symbolic},y::Unitful.AbstractQuantity) = x * y
 
 "Find the units of a symbolic item."
-get_units(x) = 1
+get_units(x::Real) = 1
 get_units(x::Unitful.Quantity) = 1 * Unitful.unit(x)
 get_units(x::Num) = get_units(value(x))
 function get_units(x::Symbolic)
@@ -53,7 +53,7 @@ function safe_get_units(term, info)
     side
 end
 
-function _validate(terms::Vector, labels::Vector; info::String = "")
+function _validate(terms::Vector, labels::Vector{String}; info::String = "")
     equnits = safe_get_units.(terms, info*", ".*labels)
     allthere = all(map(x -> x!==nothing, equnits))
     allmatching = true
@@ -86,20 +86,50 @@ function validate(eq::ModelingToolkit.Equation, noisevec::Vector; info::String =
     _validate(terms, labels, info = info)
 end
 
-function validate(eqs::Vector{ModelingToolkit.Equation})
-    all([validate(eqs[idx], info = "In eq. #$idx") for idx in 1:length(eqs)])
+function validate(eqs::Vector{ModelingToolkit.Equation}; info::String = "")
+    all([validate(eqs[idx], info = info*"in eq. #$idx") for idx in 1:length(eqs)])
 end
 
 function validate(eqs::Vector{ModelingToolkit.Equation}, noise::Vector)
-    all([validate(eqs[idx], noise[idx], info = "In eq. #$idx") for idx in 1:length(eqs)])
+    all([validate(eqs[idx], noise[idx], info = "in eq. #$idx") for idx in 1:length(eqs)])
 end
 
 function validate(eqs::Vector{ModelingToolkit.Equation}, noise::Matrix)
-    all([validate(eqs[idx], noise[idx, :], info = "In eq. #$idx") for idx in 1:length(eqs)])
+    all([validate(eqs[idx], noise[idx, :], info = "in eq. #$idx") for idx in 1:length(eqs)])
 end
 
+function validate(terms::Vector{<:SymbolicUtils.Symbolic})
+   _validate(terms,["in term #$idx" for idx in 1:length(terms)],info = "")
+end
+
+function validate(jump::Union{ModelingToolkit.VariableRateJump, ModelingToolkit.ConstantRateJump}, t::Symbolic; info::String = "")
+    _validate([jump.rate, 1/t], ["rate", "1/t"], info = info) && # Assuming the rate is per time units
+    validate(jump.affect!,info = info) 
+end
+
+function validate(jump::ModelingToolkit.MassActionJump, t::Symbolic; info::String = "")
+    left_symbols = [x[1] for x in jump.reactant_stoch] #vector of pairs of symbol,int -> vector symbols
+    net_symbols = [x[1] for x in jump.net_stoch]
+    all_symbols = vcat(left_symbols,net_symbols)
+    allgood = _validate(all_symbols, string.(all_symbols), info = info)
+    n = sum(x->x[2],jump.reactant_stoch,init = 0)
+    base_unitful = all_symbols[1] #all same, get first
+    allgood && _validate([jump.scaled_rates, 1/(t*base_unitful^n)], ["scaled_rates", "1/(t*reactants^$n))"], info = info)
+end
+
+function validate(jumps::Vector{<:JumpType}, t::Symbolic; info::String = "")
+    all([validate(jumps[idx], t, info = info*"in jump #$idx") for idx in 1:length(jumps)])
+end
+
+function validate(jumps::ArrayPartition{<:Union{Any, Vector{<:JumpType}}}, t::Symbolic)
+    labels = ["in Mass Action Jumps, ", "in Constant Rate Jumps, ", "in Variable Rate Jumps, "]
+    all([validate(jumps.x[idx], t, info = labels[idx]) for idx in 1:3])
+end
+
+validate(eqs::Vector{<:ModelingToolkit.Reaction}) = validate(oderatelaw.(eqs))
+
 "Returns true iff units of equations are valid."
-validate(eqs::Vector) = validate(convert.(ModelingToolkit.Equation, eqs))
+validate(eqs::Vector{Any}) = length(eqs) == 0 
 
 "Throws error if units of equations are invalid."
 check_units(eqs...) = validate(eqs...) || throw(ArgumentError("Some equations had invalid units. See warnings for details."))

test/units.jl

@@ -1,4 +1,4 @@
-using ModelingToolkit, Unitful, OrdinaryDiffEq
+using ModelingToolkit, Unitful, OrdinaryDiffEq, DiffEqJump
 using Test
 MT = ModelingToolkit
 @parameters τ [unit = u"ms"]
@@ -130,11 +130,11 @@ eqs = [D(E) ~ P - E/τ
 
 noiseeqs = [0.1u"MW",
             0.1u"MW"]
-sys = SDESystem(eqs,noiseeqs,t,[P,E],[τ,Q])
+sys = SDESystem(eqs, noiseeqs, t, [P, E], [τ, Q])
 # With noise matrix
 noiseeqs = [0.1u"MW" 0.1u"MW"
             0.1u"MW" 0.1u"MW"]
-sys = SDESystem(eqs,noiseeqs,t,[P,E],[τ,Q])
+sys = SDESystem(eqs,noiseeqs, t, [P, E], [τ, Q])
 
 noiseeqs = [0.1u"MW" 0.1u"MW"
             0.1u"MW" 0.1u"s"]
@@ -158,10 +158,42 @@ sys_simple = structural_simplify(sys)
 @parameters v [unit = u"m/s"] r [unit =u"m"^3/u"s"] t [unit = u"s"]
 eqs = [V ~ r*t,
        V ~ L^3]
-sys = NonlinearSystem(eqs,[V,L],[t,r])
+sys = NonlinearSystem(eqs, [V, L], [t, r])
 sys_simple = structural_simplify(sys)
 
 eqs = [L ~ v*t,
        V ~ L^3]
-sys = NonlinearSystem(eqs,[V,L],[t,r])
-sys_simple = structural_simplify(sys)
+sys = NonlinearSystem(eqs, [V,L], [t,r])
+sys_simple = structural_simplify(sys)
+
+#Jump System
+@parameters β [unit = 1/(u"mol"^2*u"s")] γ [unit = 1/(u"mol"*u"s")] t [unit = u"s"] 
+@variables S(t) [unit = u"mol"] I(t) [unit = u"mol"] R(t) [unit = u"mol"]
+rate₁   = β*S*I
+affect₁ = [S ~ S - 1u"mol", I ~ I + 1u"mol"]
+rate₂   = γ*I
+affect₂ = [I ~ I - 1u"mol", R ~ R + 1u"mol"]
+j₁      = ConstantRateJump(rate₁, affect₁)
+j₂      = VariableRateJump(rate₂, affect₂)
+js      = JumpSystem([j₁, j₂], t, [S, I, R], [β, γ])
+
+affect_wrong = [S ~ S - 1u"mol", I ~ I + 1]
+j_wrong      = ConstantRateJump(rate₁, affect_wrong)
+@test_throws ArgumentError JumpSystem([j_wrong, j₂], t, [S, I, R], [β, γ])
+
+rate_wrong   = γ^2*I
+j_wrong     = ConstantRateJump(rate_wrong, affect₂)
+@test_throws ArgumentError JumpSystem([j₁, j_wrong], t, [S, I, R], [β, γ])
+
+# mass action jump tests for SIR model
+maj1 = MassActionJump(2*β/2, [S => 1, I => 1], [S => -1, I => 1])
+maj2 = MassActionJump(γ, [I => 1], [I => -1, R => 1])
+js3  = JumpSystem([maj1, maj2], t, [S,I,R], [β,γ])
+
+#Test unusual jump system
+@parameters β γ t
+@variables S(t) I(t) R(t)
+
+maj1 = MassActionJump(2.0, [0 => 1], [S => 1])
+maj2 = MassActionJump(γ, [S => 1], [S => -1])
+js4  = JumpSystem([maj1, maj2], t, [S], [β, γ])