From 6e9bd40d1f239f7c2a987aeaae97a664d1fa608c Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 4 Apr 2024 16:57:32 -0400 Subject: [PATCH] Finish events and add test --- Project.toml | 3 +- src/reactionsystem.jl | 149 +++++++++----- test/dsl/dsl_options.jl | 89 -------- test/miscellaneous_tests/events.jl | 316 +++++++++++++++++++++++++++-- 4 files changed, 407 insertions(+), 150 deletions(-) diff --git a/Project.toml b/Project.toml index f9a66714ae..3556c93e03 100644 --- a/Project.toml +++ b/Project.toml @@ -60,6 +60,7 @@ julia = "1.9" [extras] BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665" +DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf" Graphviz_jll = "3c863552-8265-54e4-a6dc-903eb78fde85" HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" @@ -80,4 +81,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [targets] -test = ["BifurcationKit", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "NonlinearSolve", "OrdinaryDiffEq", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "StructuralIdentifiability", "Test", "Unitful"] +test = ["BifurcationKit", "DiffEqCallbacks", "DomainSets", "Graphviz_jll", "HomotopyContinuation", "NonlinearSolve", "OrdinaryDiffEq", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "StructuralIdentifiability", "Test", "Unitful"] diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index b0ef5e4904..8f66742fcd 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -628,70 +628,77 @@ function ReactionSystem(eqs, iv, unknowns, ps; continuous_events = nothing, discrete_events = nothing, metadata = nothing) - - name === nothing && + + # Error checks + if name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro")) + end sysnames = nameof.(systems) - (length(unique(sysnames)) == length(sysnames)) || - throw(ArgumentError("System names must be unique.")) + (length(unique(sysnames)) == length(sysnames)) || throw(ArgumentError("System names must be unique.")) + # Handle defaults values provided via optional arguments. if !(isempty(default_u0) && isempty(default_p)) - Base.depwarn("`default_u0` and `default_p` are deprecated. Use `defaults` instead.", - :ReactionSystem, force = true) + Base.depwarn("`default_u0` and `default_p` are deprecated. Use `defaults` instead.", :ReactionSystem, force = true) end defaults = MT.todict(defaults) defaults = Dict{Any, Any}(value(k) => value(v) for (k, v) in pairs(defaults)) + # Extracts independent variables (iv and sivs), dependent variables (species and variables) + # and parameters. Sorts so that species comes before variables in unknowns vector. iv′ = value(iv) sivs′ = if spatial_ivs === nothing Vector{typeof(iv′)}() else value.(MT.scalarize(spatial_ivs)) end - unknowns′ = sort!(value.(MT.scalarize(unknowns)), by = !isspecies) # species come first + unknowns′ = sort!(value.(MT.scalarize(unknowns)), by = !isspecies) spcs = filter(isspecies, unknowns′) ps′ = value.(MT.scalarize(ps)) + # Checks that no (by Catalyst) forbidden symbols are used. allsyms = Iterators.flatten((ps′, unknowns′)) - all(sym -> getname(sym) ∉ forbidden_symbols_error, allsyms) || + if !all(sym -> getname(sym) ∉ forbidden_symbols_error, allsyms) error("Catalyst reserves the symbols $forbidden_symbols_error for internal use. Please do not use these symbols as parameters or unknowns/species.") + end - # sort Reactions before Equations + # Handles reactions and equations. Sorts so that reactions are before equaions in the equations vector. eqs′ = CatalystEqType[eq for eq in eqs] sort!(eqs′; by = eqsortby) rxs = Reaction[rx for rx in eqs if rx isa Reaction] + # Additional error checks. if any(MT.isparameter, unknowns′) psts = filter(MT.isparameter, unknowns′) throw(ArgumentError("Found one or more parameters among the unknowns; this is not allowed. Move: $psts to be parameters.")) end - if any(isconstant, unknowns′) csts = filter(isconstant, unknowns′) throw(ArgumentError("Found one or more constant species among the unknowns; this is not allowed. Move: $csts to be parameters.")) end - - # if there are BC species, check they are balanced in their reactions + # If there are BC species, check they are balanced in their reactions. if balanced_bc_check && any(isbc, unknowns′) for rx in eqs - if rx isa Reaction - isbcbalanced(rx) || - throw(ErrorException("BC species must be balanced, appearing as a substrate and product with the same stoichiometry. Please fix reaction: $rx")) + if (rx isa Reaction) && !isbcbalanced(rx) + throw(ErrorException("BC species must be balanced, appearing as a substrate and product with the same stoichiometry. Please fix reaction: $rx")) end end end + # Adds all unknowns/parameters to the `var_to_name` vector. + # Adds their (potential) default values to the defaults vector. var_to_name = Dict() MT.process_variables!(var_to_name, defaults, unknowns′) MT.process_variables!(var_to_name, defaults, ps′) MT.collect_var_to_name!(var_to_name, eq.lhs for eq in observed) - nps = if networkproperties === nothing - NetworkProperties{Int, get_speciestype(iv′, unknowns′, systems)}() + # Computes network properties. + if networkproperties === nothing + nps = NetworkProperties{Int, get_speciestype(iv′, unknowns′, systems)}() else - networkproperties + nps = networkproperties end + # Creates the continious and discrete callbacks. ccallbacks = MT.SymbolicContinuousCallbacks(continuous_events) dcallbacks = MT.SymbolicDiscreteCallbacks(discrete_events) @@ -705,25 +712,31 @@ function ReactionSystem(rxs::Vector, iv = Catalyst.DEFAULT_IV; kwargs...) end # search the symbolic expression for parameters or unknowns -# and save in ps and sts respectively. vars is used to cache results -function findvars!(ps, sts, exprtosearch, ivs, vars) +# and save in ps and us respectively. vars is used to cache results +function findvars!(ps, us, exprtosearch, ivs, vars) MT.get_variables!(vars, exprtosearch) for var in vars (var ∈ ivs) && continue if MT.isparameter(var) push!(ps, var) else - push!(sts, var) + push!(us, var) end end empty!(vars) end -# Only used internally by the @reaction_network macro. Permits giving an initial order to -# the parameters, and then adds additional ones found in the reaction. Name could be -# changed. -function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, sts_in, ps_in; - spatial_ivs = nothing, kwargs...) +# Called internally (whether DSL-based or programmtic model creation is used). +# Creates a sorted reactions + equations vector, also ensuring reaction is first in this vector. +# Extracts potential species, variables, and parameters from the input (if not provided as part of +# the model creation) and creates the corresponding vectors. +# While species are ordered before variables in the unknowns vector, this ordering is not imposed here, +# but carried out at a later stage. +function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, us_in, ps_in; spatial_ivs = nothing, + continuous_events = [], discrete_events = [], kwargs...) + + # Creates a combined iv vector (iv and sivs). This is used later in the function (so that + # independent variables can be exluded when encountered quantities are added to `us` and `ps`). t = value(iv) ivs = Set([t]) if (spatial_ivs !== nothing) @@ -731,51 +744,93 @@ function make_ReactionSystem_internal(rxs_and_eqs::Vector, iv, sts_in, ps_in; push!(ivs, value(siv)) end end - sts = OrderedSet{eltype(sts_in)}(sts_in) + + # Initialises the new unknowns and parameter vectors. + # Preallocates the `vars` set, which is used by `findvars!` + us = OrderedSet{eltype(us_in)}(us_in) ps = OrderedSet{eltype(ps_in)}(ps_in) vars = OrderedSet() + # Extracts the reactions and equations from the combined reactions + equations input vector. all(eq -> eq isa Union{Reaction, Equation}, rxs_and_eqs) rxs = Reaction[eq for eq in rxs_and_eqs if eq isa Reaction] eqs = Equation[eq for eq in rxs_and_eqs if eq isa Equation] - # add species / parameters that are substrates / products first - for rx in rxs, reactants in (rx.substrates, rx.products) - for spec in reactants - MT.isparameter(spec) ? push!(ps, spec) : push!(sts, spec) - end - end - + # Loops through all reactions, adding encountered quantities to the unknown and parameter vectors. for rx in rxs - findvars!(ps, sts, rx.rate, ivs, vars) - for s in rx.substoich - (s isa Symbolic) && findvars!(ps, sts, s, ivs, vars) + # Loops through all reaction substrates and products, extracting these. + for reactants in (rx.substrates, rx.products), spec in reactants + MT.isparameter(spec) ? push!(ps, spec) : push!(us, spec) end - for p in rx.prodstoich - (p isa Symbolic) && findvars!(ps, sts, p, ivs, vars) + + # Adds all quantitites encountered in the reaction's rate. + findvars!(ps, us, rx.rate, ivs, vars) + + # Extracts all quantitites encountered within stoichiometries. + for stoichiometry in (rx.substoich, rx.prodstoich), sym in stoichiometry + (sym isa Symbolic) && findvars!(ps, us, sym, ivs, vars) end - end - stsv = collect(sts) - psv = collect(ps) + # Will appear here: add stuff from nosie scaling. + end + # Extracts any species, variables, and parameters that occur in (non-reaction) equations. + # Creates the new reactions + equations vector, `fulleqs` (sorted reactions first, equations next). if !isempty(eqs) osys = ODESystem(eqs, iv; name = gensym()) fulleqs = CatalystEqType[rxs; equations(osys)] - union!(stsv, unknowns(osys)) - union!(psv, parameters(osys)) + union!(us, unknowns(osys)) + union!(ps, parameters(osys)) else fulleqs = rxs - end + end - ReactionSystem(fulleqs, t, stsv, psv; spatial_ivs, kwargs...) + # Loops through all events, adding encountered quantities to the unknwon and parameter vectors. + find_event_vars!(ps, us, continuous_events, ivs, vars) + find_event_vars!(ps, us, discrete_events, ivs, vars) + + # Converts the found unknowns and parameters to vectors. + usv = collect(us) + psv = collect(ps) + + # Passes the processed input into the next `ReactionSystem` call. + ReactionSystem(fulleqs, t, usv, psv; spatial_ivs, continuous_events, discrete_events, kwargs...) end function ReactionSystem(iv; kwargs...) ReactionSystem(Reaction[], iv, [], []; kwargs...) end - +# Loops through all events in an supplied event vector, adding all unknowns and parameters found in +# its condition and affect functions to their respective vectors (`ps` and `us`). +function find_event_vars!(ps, us, events::Vector, ivs, vars) + foreach(event -> find_event_vars!(ps, us, event, ivs, vars), events) +end +# For a single event, adds quantitites from its condition and affect expression(s) to `ps` and `us`. +function find_event_vars!(ps, us, event, ivs, vars) + conds, affects = event + # For discrete events, the condition can be a single value (for periodic events). + # If not, it is a vector of conditions and we must check each. + if conds isa Vector + for cond in conds + # For continious events the conditions are equations (with lhs and rhs). + # For discrete events, they are single expressions. + if cond isa Equation + findvars!(ps, us, cond.lhs, ivs, vars) + findvars!(ps, us, cond.rhs, ivs, vars) + else + findvars!(ps, us, cond, ivs, vars) + end + end + else + findvars!(ps, us, conds, ivs, vars) + end + # The affects is always a vector of equations. Here, we handle the lhs and rhs separately. + for affect in affects + findvars!(ps, us, affect.lhs, ivs, vars) + findvars!(ps, us, affect.rhs, ivs, vars) + end +end """ remake_ReactionSystem_internal(rs::ReactionSystem; default_reaction_metadata::Vector{Pair{Symbol, T}} = Vector{Pair{Symbol, Any}}()) where {T} diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl index fed8547245..7707a9414c 100644 --- a/test/dsl/dsl_options.jl +++ b/test/dsl/dsl_options.jl @@ -860,93 +860,4 @@ let @equations X ~ p - S (P,D), 0 <--> S end -end - -### Events ### - -# Compares models with complicated events that are created programmatically/with the DSL. -# Checks that simulations are correct. -# Checks that various simulation inputs works. -# Checks continuous, discrete, preset time, and periodic events. -# Tests event affecting non-species components. - -let - # Creates model via DSL. - rn_dsl = @reaction_network rn begin - @parameters thres=1.0 dY_up - @variables Z(t) - @continuous_events begin - [t - 2.5] => [p ~ p + 0.2] - [X - thres, Y - X] => [X ~ X - 0.5, Z ~ Z + 0.1] - end - @discrete_events begin - 2.0 => [dX ~ dX + 0.1, dY ~ dY + dY_up] - [1.0, 5.0] => [p ~ p - 0.1] - (Z > Y) => [Z ~ Z - 0.1] - end - - (p, dX), 0 <--> X - (p, dY), 0 <--> Y - end - - # Creates model programmatically. - @variables t Z(t) - @species X(t) Y(t) - @parameters p dX dY thres=1.0 dY_up - rxs = [ - Reaction(p, nothing, [X], nothing, [1]) - Reaction(dX, [X], nothing, [1], nothing) - Reaction(p, nothing, [Y], nothing, [1]) - Reaction(dY, [Y], nothing, [1], nothing) - ] - continuous_events = [ - t - 2.5 => p ~ p + 0.2 - [X - thres, Y - X] => [X ~ X - 0.5, Z ~ Z + 0.1] - ] - discrete_events = [ - 2.0 => [dX ~ dX + 0.1, dY ~ dY + dY_up] - [1.0, 5.0] => [p ~ p - 0.1] - (Z > Y) => [Z ~ Z - 0.1] - ] - rn_prog = ReactionSystem(rxs, t; continuous_events, discrete_events, name=:rn) - - # Tests that approaches yield identical results. - @test isequal(rn_dsl, rn_prog) - - u0 = [X => 1.0, Y => 0.5, Z => 0.25] - tspan = (0.0, 20.0) - ps = [p => 1.0, dX => 0.5, dY => 0.5, dY_up => 0.1] - - sol_dsl = solve(ODEProblem(rn_dsl, u0, tspan, ps), Tsit5()) - sol_prog = solve(ODEProblem(rn_prog, u0, tspan, ps), Tsit5()) - @test sol_dsl == sol_prog -end - -# Compares DLS events to those given as callbacks. -# Checks that events works when given to SDEs. -let - # Creates models. - rn = @reaction_network begin - (p, d), 0 <--> X - end - rn_events = @reaction_network begin - @discrete_events begin - [5.0, 10.0] => [X ~ X + 100.0] - end - @continuous_events begin - [X ~ 90.0] => [X ~ X + 10.0] - end - (p, d), 0 <--> X - end - cb_disc = ModelingToolkit.PresetTimeCallback([5.0, 10.0], int -> (int[:X] += 100.0)) - cb_cont = ContinuousCallback((u, t, int) -> (u[1] - 90.0), int -> (int[:X] += 10.0)) - - # Simulates models. - u0 = [:X => 100.0] - tspan = (0.0, 50.0) - ps = [:p => 100.0, :d => 1.0] - sol = solve(SDEProblem(rn, u0, tspan, ps), ImplicitEM(); seed, callback = CallbackSet(cb_disc, cb_cont)) - sol_events = solve(SDEProblem(rn_events, u0, tspan, ps), ImplicitEM(); seed) - - @test sol == sol_events end \ No newline at end of file diff --git a/test/miscellaneous_tests/events.jl b/test/miscellaneous_tests/events.jl index 255f217241..df1749f47e 100644 --- a/test/miscellaneous_tests/events.jl +++ b/test/miscellaneous_tests/events.jl @@ -1,17 +1,23 @@ ### Prepares Tests ### # Fetch packages. -using Catalyst, OrdinaryDiffEq, Test +using Catalyst, DiffEqCallbacks, JumpProcesses, OrdinaryDiffEq, StochasticDiffEq, Test + +# Sets stable rng number. +using StableRNGs +rng = StableRNG(12345) +seed = rand(rng, 1:100) # Sets the default `t` and `D` to use. t = default_t() D = default_time_deriv() + ### Basic Tests ### # Test discrete event is propagated to ODE solver correctly. let - # Creates model (essentially a jagged oscillation, where `V` is reset to 1.0 every `1.0` time units). + # Creates model (essentially a jagged oscillation, where `V` is reset to 1.0 every 1.0 time units). @variables V(t)=1.0 eqs = [D(V) ~ V] discrete_events = [1.0 => [V ~ 1.0]] @@ -50,32 +56,34 @@ let osys = complete(convert(ODESystem, complete(rs))) @test length(ModelingToolkit.continuous_events(osys)) == 1 @test length(ModelingToolkit.discrete_events(osys)) == 0 - oprob = ODEProblem(rs, [], (0.0, 20.0)) + oprob = ODEProblem(osys, [], (0.0, 20.0)) sol = solve(oprob, Tsit5()) @test sol(20.0, idxs = V) ≈ 2.5 end +# Tests that species/variables/parameters only encountered in events are added to `ReactionSystem`s properly. +# Tests for both discrete and continuous events. Tests that these quantities can be accessed in Problems. +# Tests that metadata for these quantities are saved properly let # Creates model. - @parameters p d α = 1.0 - @species X(t) A(t) = 2 + @parameters p d α::Int64 = 1 + @species X(t) A(t) = 2 [description="A species"] @variables a(t) = 3 rxs = [ Reaction(p, nothing, [X]), Reaction(d, [X], nothing) ] - continuous_events = [(a == t) => A ~ A + α] - discrete_events = [2.0 => A ~ α + a] + continuous_events = [α ~ t] => [A ~ A + a] + discrete_events = [2.0 => [A ~ α + a]] @named rs_ce = ReactionSystem(rxs, t; continuous_events) @named rs_de = ReactionSystem(rxs, t; discrete_events) - continuous_events = [(a == t) => A ~ A + a] - discrete_events = [2.0 => A ~ α + a] + continuous_events = [[α ~ t] => [A ~ A + α]] + discrete_events = [2.0 => [A ~ a]] @named rs_ce_de = ReactionSystem(rxs, t; continuous_events, discrete_events) rs_ce = complete(rs_ce) rs_de = complete(rs_de) rs_ce_de = complete(rs_ce_de) - # Tests model content. issetequal(species(rs_ce), [X, A]) issetequal(species(rs_de), [X, A]) @@ -86,15 +94,297 @@ let issetequal(parameters(rs_ce), [p, d, α]) issetequal(parameters(rs_de), [p, d, α]) issetequal(parameters(rs_ce_de), [p, d, α]) + @test Symbolics.unwrap(rs_ce_de.α) isa Symbolics.BasicSymbolic{Int64} + @test Symbolics.unwrap(rs_de.α) isa Symbolics.BasicSymbolic{Int64} + @test Symbolics.unwrap(rs_ce_de.α) isa Symbolics.BasicSymbolic{Int64} + @test getdescription(rs_ce_de.A) == "A species" + @test getdescription(rs_de.A) == "A species" + @test getdescription(rs_ce_de.A) == "A species" - # Tests that species/variables/parameters can be accessed correctly one a MTK structure have been created. + # Tests that species/variables/parameters can be accessed correctly one a MTK problem have been created. u0 = [X => 1] tspan = (0.0, 10.0) ps = [p => 10.0, d => 0.2] - for XProb in [ODEProblem, SDEProblem, DiscreteProblem], rs in [rs_ce, rs_de, rs_ce_de] + for XProb in [ODEProblem, SDEProblem], rs in [rs_ce, rs_de, rs_ce_de] prob = XProb(rs, u0, (0.0, 10.0), ps) @test prob[A] == 2 @test prob[a] == 3 - @test prob.ps[α] == 1.0 + @test prob.ps[α] == 1 + @test prob.ps[α] isa Int64 + end + + # Handles `DiscreteProblem`s and `JumpProblem`s (these cannot contain continuous events or variables). + discrete_events = [2.0 => [A ~ A + α]] + @named rs_de_2 = ReactionSystem(rxs, t; discrete_events) + rs_de_2 = complete(rs_de_2) + dprob = DiscreteProblem(rs_de_2, u0, (0.0, 10.0), ps) + jprob = JumpProblem(rs_de_2, dprob, Direct()) + for prob in [dprob, jprob] + @test dprob[A] == 2 + @test dprob.ps[α] == 1 + @test dprob.ps[α] isa Int64 + end +end + + +### Event Input Checks ### + +# Checks that singular events can be provided as vectors/not as vectors and this does not matter. +let + @parameters p d + @species X(t) + rxs = [ + Reaction(p, nothing, [X]), + Reaction(d, [X], nothing) + ] + ce = [X ~ 1.0] => [X ~ 0.5] + de = [2.0] => [p ~ 1.0] + rs1 = ReactionSystem(rxs, t; continuous_events = ce, discrete_events = de, name = :rs) + rs2 = ReactionSystem(rxs, t; continuous_events = [ce], discrete_events = de, name = :rs) + rs3 = ReactionSystem(rxs, t; continuous_events = ce, discrete_events = [de], name = :rs) + rs4 = ReactionSystem(rxs, t; continuous_events = [ce], discrete_events = [de], name = :rs) + @test rs1 == rs2 == rs3 == rs4 +end + +# Checks that various various erroneous forms yield errors. +# I.e. ensures affects/conditions requires vector forms in the right cases. +let + # Prepares the model reaction. + @parameters p d + @species X(t) + rxs = [ + Reaction(p, nothing, [X]), + Reaction(d, [X], nothing) + ] + + # Declares various misformatted events . + continuous_events_bad = [ + X ~ 1.0 => [X ~ 0.5], # Scalar condition. + [X ~ 1.0] => X ~ 0.5, # Scalar affect. + (X ~ 1.0,) => [X ~ 0.5], # Tuple condition. + [X ~ 1.0] => (X ~ 0.5,), # Tuple affect. + [X - 1.0] => [X ~ 0.5], # Non-equation condition (1). + [X == 1.0] => [X ~ 0.5], # Non-equation condition (2). + # [X ~ 1.0] => [X ~ 0.5, p ~ 0.5], # No error on system creation, but permitted. Should probably throw an early error. + ] + discrete_events_bad = [ + [2.0] => p ~ 1.0, # Scalar affect. + [2.0] => (p ~ 1.0, ), # Tuple affect. + #[X > 2.0] => [p ~ 1.0], # Vector conditions. Should probably throw an error here already, currently does not. + #(1.0, 2.0) => [p ~ 1.0] # Tuple condition. Should probably throw an error here already, currently does not. + ] + + # Checks that errors are produced. + for continuous_events in continuous_events_bad + println(continuous_events) + @test_throws Exception @named rs = ReactionSystem(rxs, t; continuous_events) + end + for discrete_events in discrete_events_bad + println(discrete_events) + @test_throws Exception @named rs = ReactionSystem(rxs, t; discrete_events) + end +end + + +### DSL-based Tests ### + +# Compares models with complicated events that are created programmatically/with the DSL. +# Checks that simulations are correct. +# Checks that various simulation inputs works. +# Checks continuous, discrete, preset time, and periodic events. +# Tests event affecting non-species components. +let + # Creates model via DSL. + rn_dsl = @reaction_network rn begin + @parameters thres=1.0 dY_up + @variables Z(t) + @continuous_events begin + [t ~ 2.5] => [p ~ p + 0.2] + [X ~ thres, Y ~ X] => [X ~ X - 0.5, Z ~ Z + 0.1] + end + @discrete_events begin + 2.0 => [dX ~ dX + 0.1, dY ~ dY + dY_up] + [1.0, 5.0] => [p ~ p - 0.1] + (Z > Y) => [Z ~ Z - 0.1] + end + + (p, dX), 0 <--> X + (p, dY), 0 <--> Y + end + + # Creates model programmatically. + @variables t Z(t) + @species X(t) Y(t) + @parameters p dX dY thres=1.0 dY_up + rxs = [ + Reaction(p, nothing, [X], nothing, [1]) + Reaction(dX, [X], nothing, [1], nothing) + Reaction(p, nothing, [Y], nothing, [1]) + Reaction(dY, [Y], nothing, [1], nothing) + ] + continuous_events = [ + t - 2.5 => p ~ p + 0.2 + [X - thres, Y - X] => [X ~ X - 0.5, Z ~ Z + 0.1] + ] + discrete_events = [ + 2.0 => [dX ~ dX + 0.1, dY ~ dY + dY_up] + [1.0, 5.0] => [p ~ p - 0.1] + (Z > Y) => [Z ~ Z - 0.1] + ] + rn_prog = ReactionSystem(rxs, t; continuous_events, discrete_events, name=:rn) + + # Tests that approaches yield identical results. + @test isequal(rn_dsl, rn_prog) + + u0 = [X => 1.0, Y => 0.5, Z => 0.25] + tspan = (0.0, 20.0) + ps = [p => 1.0, dX => 0.5, dY => 0.5, dY_up => 0.1] + + sol_dsl = solve(ODEProblem(rn_dsl, u0, tspan, ps), Tsit5()) + sol_prog = solve(ODEProblem(rn_prog, u0, tspan, ps), Tsit5()) + @test sol_dsl == sol_prog + + sol_dsl = solve(SDEProblem(rn_dsl, u0, tspan, ps), ImplicitEM(); seed = 1234) + sol_prog = solve(SDEProblem(rn_prog, u0, tspan, ps), ImplicitEM(); seed = 1234) + @test sol_dsl == sol_prog +end + +# Checks that misformatted events yields errors in the DSL. +let + # Quantity in event not declared elsewhere (continuous events). + @test_throws Exception @eval @reaction_network begin + @continuous_events X ~ 2.0 => [X ~ X + 1] + end + + # Scalar condition (continuous events). + @test_throws Exception @eval @reaction_network begin + @species X(t) + @continuous_events X ~ 2.0 => [X ~ X + 1] + end + + # Scalar affect (continuous events). + @test_throws Exception @eval @reaction_network begin + @species X(t) + @continuous_events [X ~ 2.0] => X ~ X + 1 + end + + # Tuple condition (continuous events). + @test_throws Exception @eval @reaction_network begin + @species X(t) + @continuous_events (X ~ 2.0,) => [X ~ X + 1] + end + + # Tuple affect (continuous events). + @test_throws Exception @eval @reaction_network begin + @species X(t) + @continuous_events [X ~ 2.0] => (X ~ X + 1,) + end + + # Non-equation condition (continuous events). + @test_throws Exception @eval @reaction_network begin + @species X(t) + @continuous_events [X - 2.0] => [X ~ X + 1] + end + + # Quantity in event not declared elsewhere (discrete events). + @test_throws Exception @eval @reaction_network begin + @discrete_events X ~ 2.0 => [X ~ X + 1] + end + + # Scalar affect (discrete events). + @test_throws Exception @eval @reaction_network begin + @species X(t) + @discrete_events 1.0 => X ~ X + 1 + end + + # Tuple affect (discrete events). + @test_throws Exception @eval @reaction_network begin + @species X(t) + @discrete_events 1.0 => (X ~ X + 1, ) + end + + # Equation condition (discrete events). + @test_throws Exception @eval @reaction_network begin + @species X(t) + @discrete_events X ~ 1.0 => [X ~ X + 1] + end +end + + +### Additional Correctness Tests ### + +# Compares simulations using MTK type events with those generated through callbacks. +# Jump simulations must be handles differently (since these only accepts discrete callbacks). +# Checks for all types of discrete callbacks, and for continuous callbacks. +# Turns of noise for SDE simulations (not sure seeding works when callbacks/events declared differently). +let + # Creates models. Jump simulations requires one with discrete events only. + rn = @reaction_network begin + @default_noise_scaling 0.0 + @parameters add::Int64 + (p,d), 0 <--> X + (p,d), 0 <--> Y + end + rn_events = @reaction_network begin + @default_noise_scaling 0.0 + @parameters add::Int64 + @continuous_events begin + [X ~ 90.0] => [X ~ X + 10.0] + end + @discrete_events begin + [5.0, 10.0] => [X ~ X + add, Y ~ Y + add] + 20.0 => [X ~ X + add] + (Y < X) => [Y ~ Y + add] + end + (p,d), 0 <--> X + (p,d), 0 <--> Y + end + rn_dics_events = @reaction_network begin + @parameters add::Int64 + @discrete_events begin + [5.0, 10.0] => [X ~ X + add, Y ~ Y + add] + 20.0 => [X ~ X + add] + (Y < X) => [Y ~ Y + add] + end + (p,d), 0 <--> X + (p,d), 0 <--> Y + end + + # Sets simulation inputs. + u0 = [:X => 100, :Y => 150] + tspan = (0.0, 90.0) + ps = [:p => 100.0, :d => 1.0, :add => 100] + + # Create callbacks + cb_cont = ContinuousCallback((u, t, int) -> (int[:X] - 90.0), int -> (int[:X] += 10.0)) + cb_disc_1 = PresetTimeCallback([5.0, 10.0], int -> (int[:X] += int.ps[:add]; int[:Y] += int.ps[:add];)) + cb_disc_2 = PresetTimeCallback(20.0:20.0:tspan[end], int -> (int[:X] += int.ps[:add])) + cb_disc_3 = DiscreteCallback((u,t,i) -> i[:Y] < i[:X], int -> (int[:Y] += int.ps[:add])) + callback = CallbackSet(cb_cont, cb_disc_1, cb_disc_2, cb_disc_3) + + # Checks for ODE simulations. + oprob = ODEProblem(rn, u0, tspan, ps) + oprob_events = ODEProblem(rn_events, u0, tspan, ps) + osol = solve(oprob, Tsit5(); callback) + osol_events = solve(oprob_events, Tsit5()) + @test osol == osol_events + + # Checks for SDE simulations. + sprob = SDEProblem(rn, u0, tspan, ps) + sprob_events = SDEProblem(rn_events, u0, tspan, ps) + ssol = solve(sprob, ImplicitEM(); seed, callback) + ssol_events = solve(sprob_events, ImplicitEM(); seed) + @test ssol == ssol_events + + # Checks for Jump simulations. + callback = CallbackSet(cb_disc_1, cb_disc_2, cb_disc_3) + dprob = DiscreteProblem(rn, u0, tspan, ps) + dprob_events = DiscreteProblem(rn_dics_events, u0, tspan, ps) + jprob = JumpProblem(rn, dprob, Direct(); rng) + jprob_events = JumpProblem(rn_dics_events, dprob_events, Direct(); rng) + sol = solve(jprob, SSAStepper(); seed, callback) + @test_broken let # Broken due to. Even if fixed, seeding might not work due to events. + sol_events = solve(jprob_events, SSAStepper(); seed) + @test sol == sol_events end end \ No newline at end of file