add ReactionSystem and start SDESystem components

Author: ChrisRackauckas

src/ModelingToolkit.jl

@@ -16,6 +16,7 @@ using MacroTools
 import MacroTools: splitdef, combinedef, postwalk, striplines
 import GeneralizedGenerated
 using DocStringExtensions
+using Base: RefValue
 
 """
 $(TYPEDEF)
@@ -91,7 +92,7 @@ include("simplify.jl")
 include("utils.jl")
 include("direct.jl")
 include("domains.jl")
-include("systems/diffeqs/diffeqsystem.jl")
+include("systems/diffeqs/odesystem.jl")
 include("systems/diffeqs/first_order_transform.jl")
 include("systems/nonlinear/nonlinear_system.jl")
 include("systems/pde/pdesystem.jl")

src/latexify_recipes.jl

@@ -22,7 +22,7 @@
     return lhs, rhs
 end
 
-@latexrecipe function f(eqs::Vector{ModelingToolkit.DiffEq}; iv=:t)
+@latexrecipe function f(eqs::Vector{ModelingToolkit.ODEExpr}; iv=:t)
     # Set default option values.
     env --> :align
 

src/systems/diffeqs/diagonalsdesystem.jl

@@ -1,9 +1,5 @@
 export ODESystem, ODEFunction
 
-
-using Base: RefValue
-
-
 isintermediate(eq::Equation) = !(isa(eq.lhs, Operation) && isa(eq.lhs.op, Differential))
 
 function flatten_differential(O::Operation)
@@ -14,21 +10,25 @@ function flatten_differential(O::Operation)
     return (x, t, order + 1)
 end
 
+<<<<<<< master:src/systems/diffeqs/diffeqsystem.jl
 struct DiffEq  # dⁿx/dtⁿ = rhs
+=======
+struct ODEExpr  # dⁿx/dtⁿ = rhs
+>>>>>>> add ReactionSystem and start SDESystem components:src/systems/diffeqs/odesystem.jl
     x::Variable
     n::Int
     rhs::Expression
 end
-function to_diffeq(eq::Equation)
+function Base.convert(::Type{ODEExpr},eq::Equation)
     isintermediate(eq) && throw(ArgumentError("intermediate equation received"))
     (x, t, n) = flatten_differential(eq.lhs)
     (isa(t, Operation) && isa(t.op, Variable) && isempty(t.args)) ||
         throw(ArgumentError("invalid independent variable $t"))
     (isa(x, Operation) && isa(x.op, Variable) && length(x.args) == 1 && isequal(first(x.args), t)) ||
         throw(ArgumentError("invalid dependent variable $x"))
-    return t.op, DiffEq(x.op, n, eq.rhs)
+    return t.op, ODEExpr(x.op, n, eq.rhs)
 end
-Base.:(==)(a::DiffEq, b::DiffEq) = isequal((a.x, a.n, a.rhs), (b.x, b.n, b.rhs))
+Base.:(==)(a::ODEExpr, b::ODEExpr) = isequal((a.x, a.n, a.rhs), (b.x, b.n, b.rhs))
 
 """
 $(TYPEDEF)
@@ -56,7 +56,7 @@ de = ODESystem(eqs)
 """
 struct ODESystem <: AbstractSystem
     """The ODEs defining the system."""
-    eqs::Vector{DiffEq}
+    eqs::Vector{ODEExpr}
     """Independent variable."""
     iv::Variable
     """Dependent (state) variables."""
@@ -86,7 +86,7 @@ struct ODESystem <: AbstractSystem
 end
 
 function ODESystem(eqs)
-    reformatted = to_diffeq.(eqs)
+    reformatted = convert.(ODEExpr,eqs)
 
     ivs = unique(r[1] for r ∈ reformatted)
     length(ivs) == 1 || throw(ArgumentError("one independent variable currently supported"))
@@ -102,8 +102,12 @@ function ODESystem(eqs)
     ODESystem(deqs, iv, dvs, ps)
 end
 
+<<<<<<< master:src/systems/diffeqs/diffeqsystem.jl
 function ODESystem(deqs::AbstractVector{DiffEq}, iv, dvs, ps)
     tgrad = RefValue(Vector{Expression}(undef, 0))
+=======
+function ODESystem(deqs::AbstractVector{ODEExpr}, iv, dvs, ps)
+>>>>>>> add ReactionSystem and start SDESystem components:src/systems/diffeqs/odesystem.jl
     jac = RefValue(Matrix{Expression}(undef, 0, 0))
     Wfact   = RefValue(Matrix{Expression}(undef, 0, 0))
     Wfact_t = RefValue(Matrix{Expression}(undef, 0, 0))
@@ -114,7 +118,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps)
     _dvs = [deq.op for deq ∈ dvs]
     _iv = iv.op
     _ps = [p.op for p ∈ ps]
-    ODESystem(getindex.(to_diffeq.(deqs),2), _iv, _dvs, _ps)
+    ODESystem(getindex.(convert.(ODEExpr,deqs),2), _iv, _dvs, _ps)
 end
 
 function _eq_unordered(a, b)

src/systems/diffeqs/diffeqsystem.jl src/systems/diffeqs/odesystem.jlsrc/systems/diffeqs/diffeqsystem.jl renamed to src/systems/diffeqs/odesystem.jl

@@ -0,0 +1,31 @@
+struct Reaction
+    rate
+    reactants
+    products
+end
+
+struct ReactionSystem
+    rxs
+    iv
+    dvs
+    ps
+end
+
+function Base.convert(::Type{<:ODESystem},rs::ReactionSystem)
+    D = Differential(rs.iv)
+    eqs = [D(x) ~ 0 for x in rs.dvs]
+
+    for rx in rs.rxs
+        for reactant in rx.reactants
+            i = findfirst(x->x.op == reactant.op,rs.dvs)
+            eqs[i] = Equation(eqs[i].lhs,eqs[i].rhs - reactant)
+        end
+
+        for product in rx.products
+            i = findfirst(x->x.op == product.op,rs.dvs)
+            eqs[i] = Equation(eqs[i].lhs,+(eqs[i].rhs,rx.reactants...))
+        end
+    end
+
+    ODESystem(eqs,rs.iv,rs.dvs,rs.ps)
+end

src/systems/reaction/reactionsystem.jl

@@ -0,0 +1,8 @@
+using ModelingToolkit
+
+@parameters t a b
+@variables x(t) y(t) z(t)
+rxs = [Reaction(a,[x],[y]),
+       Reaction(b,[y],[x])]
+rs = ReactionSystem(rxs,t,[x,y],[a,b])
+sys = convert(ODESystem,rs)