Merge pull request #1187 from isaacsas/move-reaction-sys-to-catalyst

mv ReactionSystems to Catalyst and fix dependencies for JumpSystems

Author: isaacsas

docs/make.jl

@@ -40,7 +40,6 @@ makedocs(
             "systems/NonlinearSystem.md",
             "systems/OptimizationSystem.md",
             "systems/ControlSystem.md",
-            "systems/ReactionSystem.md",
             "systems/PDESystem.md",
         ],
         "comparison.md",

docs/src/basics/AbstractSystem.md

@@ -4,7 +4,7 @@
 
 The `AbstractSystem` interface is the core of the system level of ModelingToolkit.jl.
 It establishes a common set of functionality that is used between systems
-from ODEs and chemical reactions, allowing users to have a common framework for
+representing ODEs, PDEs, SDEs and more, allowing users to have a common framework for
 model manipulation and compilation.
 
 ### Subtypes

docs/src/index.md

@@ -72,7 +72,7 @@ is built on, consult the
 - Nonlinear systems
 - Optimization problems
 - Continuous-Time Markov Chains
-- Chemical Reactions
+- Chemical Reactions (via [Catalyst.jl](https://github.com/SciML/Catalyst.jl))
 - Nonlinear Optimal Control
 
 ## Model Import Formats

docs/src/systems/ReactionSystem.md

@@ -135,7 +135,6 @@ 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")
@@ -172,7 +171,6 @@ export alias_elimination, flatten, connect, @connector
 export ode_order_lowering, liouville_transform
 export runge_kutta_discretize
 export PDESystem
-export Reaction, ReactionSystem, ismassaction, oderatelaw, jumpratelaw
 export Differential, expand_derivatives, @derivatives
 export Equation, ConstrainedEquation
 export Term, Sym

src/ModelingToolkit.jl

@@ -13,32 +13,25 @@ Notes:
 Example:
 ```julia
 using ModelingToolkit
-@parameters β γ κ η t
-@variables S(t) I(t) R(t)
+@parameters β γ κ η 
+@variables t S(t) I(t) R(t)
 
-# use a reaction system to easily generate ODE and jump systems
-rxs = [Reaction(β, [S,I], [I], [1,1], [2]),
-       Reaction(γ, [I], [R]),
-       Reaction(κ+η, [R], [S])]
-rs = ReactionSystem(rxs, t, [S,I,R], [β,γ,κ,η])
+rate₁   = β*S*I
+rate₂   = γ*I+t
+affect₁ = [S ~ S - 1, I ~ I + 1]
+affect₂ = [I ~ I - 1, R ~ R + 1]
+j₁ = ConstantRateJump(rate₁,affect₁)
+j₂ = VariableRateJump(rate₂,affect₂)
 
-# ODEs:
-odesys = convert(ODESystem, rs)
-
-# dependency of each ODE on state variables
-equation_dependencies(odesys)
-
-# dependency of each ODE on parameters
-equation_dependencies(odesys, variables=parameters(odesys))
-
-# Jumps
-jumpsys = convert(JumpSystem, rs)
+# create a JumpSystem using these jumps
+@named jumpsys = JumpSystem([j₁,j₂], t, [S,I,R], [β,γ])
 
 # dependency of each jump rate function on state variables
 equation_dependencies(jumpsys)
 
 # dependency of each jump rate function on parameters
 equation_dependencies(jumpsys, variables=parameters(jumpsys))
+
 ```
 """
 function equation_dependencies(sys::AbstractSystem; variables=states(sys))

src/systems/dependency_graphs.jl

@@ -165,11 +165,12 @@ end
 
 function numericrstoich(mtrs::Vector{Pair{V,W}}, statetoid) where {V,W}
     rs = Vector{Pair{Int,W}}()
-    for (spec,stoich) in mtrs
+    for (wspec,stoich) in mtrs
+        spec = value(wspec)
         if !istree(spec) && _iszero(spec)
             push!(rs, 0 => stoich)
         else
-            push!(rs, statetoid[value(spec)] => stoich)
+            push!(rs, statetoid[spec] => stoich)
         end
     end
     sort!(rs)
@@ -178,7 +179,8 @@ end
 
 function numericnstoich(mtrs::Vector{Pair{V,W}}, statetoid) where {V,W}
     ns = Vector{Pair{Int,W}}()
-    for (spec,stoich) in mtrs
+    for (wspec,stoich) in mtrs
+        spec = value(wspec)
         !istree(spec) && _iszero(spec) && error("Net stoichiometry can not have a species labelled 0.")
         push!(ns, statetoid[spec] => stoich)
     end
@@ -301,7 +303,8 @@ end
 
 ### Functions to determine which states a jump depends on
 function get_variables!(dep, jump::Union{ConstantRateJump,VariableRateJump}, variables)
-    (jump.rate isa Symbolic) && get_variables!(dep, jump.rate, variables)
+    jr = value(jump.rate)
+    (jr isa Symbolic) && get_variables!(dep, jr, variables)
     dep
 end