save progress

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

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

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 ###############################
 
 """

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
 

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...)

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,15 +739,15 @@ 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
         end 
         (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

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