From f362a3724bed0639e33a17256155693feb803451 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 2 Apr 2024 20:32:10 -0400 Subject: [PATCH] save progress --- docs/src/api/catalyst_api.md | 11 --- src/Catalyst.jl | 5 +- src/networkapi.jl | 87 ------------------- src/reaction_network.jl | 6 +- src/reactionsystem.jl | 13 ++- test/dsl/dsl_options.jl | 29 ++++--- test/miscellaneous_tests/events.jl | 77 +++++++++++++--- .../hybrid_equation_reaction_systems.jl | 0 8 files changed, 97 insertions(+), 131 deletions(-) create mode 100644 test/reactionsystem_structure/hybrid_equation_reaction_systems.jl diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md index 1cfd8abf60..5085a60d3f 100644 --- a/docs/src/api/catalyst_api.md +++ b/docs/src/api/catalyst_api.md @@ -157,17 +157,6 @@ species nonspecies reactionparams reactions -is_reaction -is_alg_equation -is_diff_equation -alg_equations -diff_equations -has_alg_equations -has_diff_equations -get_alg_eqs -get_diff_eqs -has_alg_eqs -has_diff_eqs numspecies numparams numreactions diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 5fc467d4f2..3e82b2cf73 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -35,11 +35,15 @@ import ModelingToolkit: check_variables, check_parameters, _iszero, _merge, check_units, get_unit, check_equations, iscomplete +# Event-related MTK stuff. +import ModelingToolkit: PresetTimeCallback + import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show import MacroTools, Graphs import Graphs: DiGraph, SimpleGraph, SimpleDiGraph, vertices, edges, add_vertices!, nv, ne import DataStructures: OrderedDict, OrderedSet import Parameters: @with_kw_noshow +import Symbolics: occursin, wrap # globals for the modulate function default_time_deriv() @@ -86,7 +90,6 @@ export mm, mmr, hill, hillr, hillar # functions to query network properties include("networkapi.jl") export species, nonspecies, reactionparams, reactions, speciesmap, paramsmap -export has_diff_equations, diff_equations, has_alg_equations, alg_equations export numspecies, numreactions, numreactionparams, setdefaults!, symmap_to_varmap export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsmap export dependants, dependents, substoichmat, prodstoichmat, netstoichmat diff --git a/src/networkapi.jl b/src/networkapi.jl index 72f385d46d..0cd3c2bdf8 100644 --- a/src/networkapi.jl +++ b/src/networkapi.jl @@ -586,93 +586,6 @@ symmap_to_varmap(sys, symmap) = symmap #error("symmap_to_varmap requires a Dict, AbstractArray or Tuple to map Symbols to values.") -### Equation Handling Accessors ### - -""" - is_reaction(eq::Equation) - -Returns `true` if the input is a `Reaction`, else false. -""" -function is_reaction(eq::Equation) - return eq isa Reaction -end - -""" - is_alg_equation(eq::Equation) - -Returns `true` if the input is an algebraic equation that does not contain any differentials. -""" -function is_alg_equation(eq) - return isdefined(eq, :lhs) && isdefined(eq, :rhs) && !is_diff_equation(eq) -end - -""" - is_diff_equation(eq::Equation) - -Returns `true` if the input is an that contains at least one differential. -""" -function is_diff_equation(eq) - isdefined(eq, :lhs) && occursin(is_derivative, wrap(eq.lhs)) && (return true) - isdefined(eq, :rhs) && occursin(is_derivative, wrap(eq.rhs)) && (return true) - return false -end - -""" - alg_equations(rs) - -Returns all the algebraic equations (that does not contain differnetials) in `rs` and its subsystems. -""" -alg_equations(rs) = filter(is_alg_equation, equations(rs)) - -""" - diff_equations(rs) - -Returns all the differnetials equations (equations that contain differnetials) in `rs` and its subsystems. -""" -diff_equations(rs) = filter(is_diff_equation, equations(rs)) - -""" - has_alg_equations(rs) - -Returns true if `rs` (or any of its subsystems) has an algebraic equation (that does not contain a differential). -""" -has_alg_equations(rs) = any(is_alg_equation, equations(rs)) - -""" - has_diff_equations(rs) - -Returns true if `rs` (or any of its subsystems) has a differnetials equations (equations that contain differnetials) . -""" -has_diff_equations(rs) = any(is_diff_equation, equations(rs)) - -""" - get_alg_eqs(rs) - -Returns all the algebraic equations (that does not contain differnetials) in `rs` and (but not its subsystems). -""" -get_alg_eqs(rs) = filter(is_alg_equation, get_eqs(rs)) - -""" - diff_equations(rs) - -Returns all the differnetial equations (equations that contain differnetials) in `rs` (but not its subsystems). -""" -get_diff_eqs(rs) = filter(is_diff_equation, get_eqs(rs)) - -""" - has_alg_equations(rs) - -Returns true if `rs` has an algebraic equation (that does not contain a differential). Does not consider subsystems. -""" -has_alg_eqs(rs) = any(is_alg_equation, get_eqs(rs)) - -""" - has_diff_equations(rs) - -Returns true if `rs` has a differnetial equations (equations that contain differnetials). Does not consider subsystems. -""" -has_diff_eqs(rs) = any(is_diff_equation, get_eqs(rs)) - ######################## reaction complexes and reaction rates ############################### """ diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 23ddb7ea2a..57f6fd1ed8 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -879,13 +879,13 @@ function read_events_option(options, event_type::Symbol) for arg in events_input.args # Formatting error checks. # NOTE: Maybe we should move these deeper into the system (rather than the DSL), throwing errors more generally? - if (arg.head != :call) || (arg.args[1] != :(=>)) || length(arg.args) != 3 + if (arg isa Expr) && (arg.head != :call) || (arg.args[1] != :(=>)) || length(arg.args) != 3 error("Events should be on form `condition => affect`, separated by a `=>`. This appears not to be the case for: $(arg).") end - if (arg.args[2] != :call) && (event_type == :continuous_events) + if (arg isa Expr) && (arg.args[2] isa Expr) && (arg.args[2].head != :vect) && (event_type == :continuous_events) error("The condition part of continious events (the left-hand side) must be a vector. This is not the case for: $(arg).") end - if arg.args[3] != :call + if (arg isa Expr) && (arg.args[3] isa Expr) && (arg.args[3].head != :vect) error("The affect part of all events (the righ-hand side) must be a vector. This is not the case for: $(arg).") end diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 088eb0cde0..fac41cad30 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -309,6 +309,12 @@ function isbcbalanced(rx::Reaction) true end +### Reaction Acessor Functions ### + +# Overwrites functions in ModelingToolkit to give the correct input. +is_diff_equation(rx::Reaction) = false +is_alg_equation(rx::Reaction) = false + ################################## Reaction Complexes #################################### """ @@ -1672,13 +1678,14 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan, pmap = symmap_to_varmap(rs, p) osys = convert(ODESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks, remove_conserved) - osys = complete(osys) # Handles potential Differential algebraic equations. if structural_simplify (osys = MT.structural_simplify(osys)) - # elseif has_alg_equations(rs) - # error("The input ReactionSystem has algebraic equations. This requires setting `structural_simplify=true` within `ODEProblem` call.") + elseif has_alg_equations(rs) + error("The input ReactionSystem has algebraic equations. This requires setting `structural_simplify=true` within `ODEProblem` call.") + else + osys = complete(osys) end return ODEProblem(osys, u0map, tspan, pmap, args...; check_length, kwargs...) diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl index 2c0e7b0dbe..e8e6c48487 100644 --- a/test/dsl/dsl_options.jl +++ b/test/dsl/dsl_options.jl @@ -3,7 +3,12 @@ ### Prepares Tests ### # Fetch packages. -using Catalyst, ModelingToolkit, OrdinaryDiffEq, Plots, Test +using Catalyst, ModelingToolkit, OrdinaryDiffEq, StochasticDiffEq, Plots, Test + +# Sets stable rng number. +using StableRNGs +rng = StableRNG(12345) +seed = rand(rng, 1:100) # Sets the default `t` to use. t = default_t() @@ -725,7 +730,7 @@ end # Checks hierarchical model. let - base_rn = @reaction_network begin + base_rn = @network_component begin @variables V1(t) @equations begin X*3V1 ~ X - 2 @@ -734,7 +739,7 @@ let end @unpack X, V1, p, d = base_rn - internal_rn = @reaction_network begin + internal_rn = @network_component begin @variables V2(t) @equations begin X*4V2 ~ X - 3 @@ -742,7 +747,7 @@ let (p,d), 0 <--> X end - rn = compose(base_rn, [internal_rn]) + rn = complete(compose(base_rn, [internal_rn])) u0 = [V1 => 1.0, X => 3.0, internal_rn.V2 => 2.0, internal_rn.X => 4.0] ps = [p => 1.0, d => 0.2, internal_rn.p => 2.0, internal_rn.d => 0.5] @@ -832,7 +837,7 @@ let @discrete_events begin 2.0 => [dX ~ dX + 0.1, dY ~ dY + dY_up] [1.0, 5.0] => [p ~ p - 0.1] - [Z > Y, Z > X] => [Z ~ Z - 0.1] + (Z > Y) => [Z ~ Z - 0.1] end (p, dX), 0 <--> X @@ -856,9 +861,9 @@ let discrete_events = [ 2.0 => [dX ~ dX + 0.1, dY ~ dY + dY_up] [1.0, 5.0] => [p ~ p - 0.1] - [Z > Y, Z > X] => [Z ~ Z - 0.1] + (Z > Y) => [Z ~ Z - 0.1] ] - rn_prog = ReactionSystem([rx1, rx2, eq], t; continuous_events, discrete_events, name=:rn) + rn_prog = ReactionSystem(rxs, t; continuous_events, discrete_events, name=:rn) # Tests that approaches yield identical results. @test isequal(rn_dsl, rn_prog) @@ -884,19 +889,19 @@ let [5.0, 10.0] => [X ~ X + 100.0] end @continuous_events begin - [X - 90.0] => [X ~ X + 10.0] + [X ~ 90.0] => [X ~ X + 10.0] end (p, d), 0 <--> X end - cb_disc = PresetTimeCallback([5.0, 10.0], int -> (int[:X] += 100.0)) + 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,. + # 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 = 1234, callback = CallbackSet(cb_disc, cb_cont)) - sol_events = solve(SDEProblem(rn_events, u0, tspan, ps), ImplicitEM(); seed = 1234) + 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 76e51d3fc3..255f217241 100644 --- a/test/miscellaneous_tests/events.jl +++ b/test/miscellaneous_tests/events.jl @@ -11,41 +11,90 @@ D = default_time_deriv() # Test discrete event is propagated to ODE solver correctly. let - # Creates model. + # 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]] - rxs = [(@reaction $V, 0 --> A), (@reaction 1.0, A --> 0)] + rxs = [ + @reaction $V, 0 --> A + @reaction 1.0, A --> 0 + ] @named rs = ReactionSystem([rxs; eqs], t; discrete_events) - rs = complete(rs) - @test length(ModelingToolkit.discrete_events(rs)) == 1 @test length(ModelingToolkit.continuous_events(rs)) == 0 + @test length(ModelingToolkit.discrete_events(rs)) == 1 # Tests in simulation. - setdefaults!(rs, [:A => 0.0]) - osys = complete(convert(ODESystem, rs)) + osys = complete(convert(ODESystem, complete(rs))) + @test length(ModelingToolkit.continuous_events(osys)) == 0 @test length(ModelingToolkit.discrete_events(osys)) == 1 - oprob = ODEProblem(osys, [], (0.0, 20.0)) + oprob = ODEProblem(osys, [osys.A => 0.0], (0.0, 20.0)) sol = solve(oprob, Tsit5()) - @test sol(10 + 10 * eps(), idxs = V) ≈ 1.0 + @test sol(10 + 10*eps(), idxs = V) ≈ 1.0 end # Test continuous event is propagated to the ODE solver. let - # Creates model. + # Creates model (a production/degradation system, but both reactions stop at `t=2.5`). @parameters α=5.0 β=1.0 - @species V(t) = 0.0 - rxs = [Reaction(α, nothing, [V]), Reaction(β, [V], nothing)] + @species V(t)=0.0 + rxs = [ + Reaction(α, nothing, [V]), + Reaction(β, [V], nothing) + ] continuous_events = [V ~ 2.5] => [α ~ 0, β ~ 0] @named rs = ReactionSystem(rxs, t; continuous_events) - rs = complete(rs) - @test length(ModelingToolkit.discrete_events(rs)) == 0 @test length(ModelingToolkit.continuous_events(rs)) == 1 + @test length(ModelingToolkit.discrete_events(rs)) == 0 # Tests in simulation. - osys = convert(ODESystem, rs) + 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)) sol = solve(oprob, Tsit5()) @test sol(20.0, idxs = V) ≈ 2.5 end + +let + # Creates model. + @parameters p d α = 1.0 + @species X(t) A(t) = 2 + @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] + @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] + @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]) + issetequal(species(rs_ce_de), [X, A]) + issetequal(unknowns(rs_ce), [X, A, a]) + issetequal(unknowns(rs_de), [X, A, a]) + issetequal(unknowns(rs_ce_de), [X, A, a]) + issetequal(parameters(rs_ce), [p, d, α]) + issetequal(parameters(rs_de), [p, d, α]) + issetequal(parameters(rs_ce_de), [p, d, α]) + + # Tests that species/variables/parameters can be accessed correctly one a MTK structure 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] + prob = XProb(rs, u0, (0.0, 10.0), ps) + @test prob[A] == 2 + @test prob[a] == 3 + @test prob.ps[α] == 1.0 + end +end \ No newline at end of file diff --git a/test/reactionsystem_structure/hybrid_equation_reaction_systems.jl b/test/reactionsystem_structure/hybrid_equation_reaction_systems.jl new file mode 100644 index 0000000000..e69de29bb2