Merge pull request #284 from SciML/reaction

Setup ReactionSystem on AbstractSystem interface

Author: ChrisRackauckas

README.md

@@ -296,15 +296,15 @@ and a parameters `β₁` and `β₂`.
 
 
 ```julia
-t = Variable(:t; known = true)()  # independent variables are treated as known
-α = Variable(:α; known = true)()  # parameters are known
-σ = Variable(:σ; known = true)    # left uncalled, since it is used as a function
-w = Variable(:w; known = false)   # unknown, left uncalled
-x = Variable(:x; known = false)(t)  # unknown, depends on `t`
-y = Variable(:y; known = false)()   # unknown, no dependents
-z = Variable(:z; known = false)(t, α, x)  # unknown, multiple arguments
-β₁ = Variable(:β, 1; known = true)() # with index 1
-β₂ = Variable(:β, 2; known = true)() # with index 2
+t = Variable(:t)()  # independent variables are treated as known
+α = Variable(:α)()  # parameters are known
+σ = Variable(:σ)    # left uncalled, since it is used as a function
+w = Variable(:w)   # unknown, left uncalled
+x = Variable(:x)(t)  # unknown, depends on `t`
+y = Variable(:y)()   # unknown, no dependents
+z = Variable(:z)(t, α, x)  # unknown, multiple arguments
+β₁ = Variable(:β, 1)() # with index 1
+β₂ = Variable(:β, 2)() # with index 2
 
 expr = β₁ * x + y^α + σ(3) * (z - t) - β₂ * w(t - 1)
 ```

src/systems/reaction/reactionsystem.jl

@@ -1,48 +1,56 @@
 struct Reaction
     rate
-    reactants
-    products
+    reactants::Vector{Operation}
+    products::Vector{Operation}
 end
 
-struct ReactionSystem
-    rxs
-    iv
-    dvs
-    ps
+struct ReactionSystem <: AbstractSystem
+    eqs::Vector{Reaction}
+    iv::Variable
+    states::Vector{Variable}
+    ps::Vector{Variable}
+    name::Symbol
+    systems::Vector{ReactionSystem}
+end
+
+function ReactionSystem(eqs,iv,dvs,ps;
+                        systems = ReactionSystem[],
+                        name=gensym(:ReactionSystem))
+    ReactionSystem(eqs,iv,convert.(Variable,dvs),convert.(Variable,ps),name,systems)
 end
 
 # TODO: Make it do the combinatorics stuff
 reaction_expr(reactants) = *(reactants...)
 
 function essemble_drift(rs)
-    D = Differential(rs.iv)
-    eqs = [D(x) ~ 0 for x in rs.dvs]
+    D = Differential(rs.iv())
+    eqs = [D(x(rs.iv())) ~ 0 for x in rs.states]
 
-    for rx in rs.rxs
+    for rx in rs.eqs
         for reactant in rx.reactants
-            i = findfirst(x->x.op == reactant.op,rs.dvs)
+            i = findfirst(x->x == reactant.op,rs.states)
             eqs[i] = Equation(eqs[i].lhs,eqs[i].rhs - rx.rate * reaction_expr(rx.reactants))
         end
 
         for product in rx.products
-            i = findfirst(x->x.op == product.op,rs.dvs)
+            i = findfirst(x->x == product.op,rs.states)
             eqs[i] = Equation(eqs[i].lhs,eqs[i].rhs + rx.rate * reaction_expr(rx.reactants))
         end
     end
     eqs
 end
 
 function essemble_diffusion(rs)
-    eqs = Expression[Constant(0) for x in rs.dvs, y in rs.rxs]
+    eqs = Expression[Constant(0) for x in rs.states, y in rs.eqs]
 
-    for (j,rx) in enumerate(rs.rxs)
+    for (j,rx) in enumerate(rs.eqs)
         for reactant in rx.reactants
-            i = findfirst(x->x.op == reactant.op,rs.dvs)
+            i = findfirst(x->x == reactant.op,rs.states)
             eqs[i,j] = -sqrt(rx.rate) * reaction_expr(rx.reactants)
         end
 
         for product in rx.products
-            i = findfirst(x->x.op == product.op,rs.dvs)
+            i = findfirst(x->x == product.op,rs.states)
             eqs[i,j] = sqrt(rx.rate) * reaction_expr(rx.reactants)
         end
     end
@@ -51,12 +59,13 @@ end
 
 function Base.convert(::Type{<:ODESystem},rs::ReactionSystem)
     eqs = essemble_drift(rs)
-    ODESystem(eqs,rs.iv,rs.dvs,rs.ps)
+    ODESystem(eqs,rs.iv,rs.states,rs.ps,name=rs.name,
+              systems=convert.(ODESystem,rs.systems))
 end
 
 function Base.convert(::Type{<:SDESystem},rs::ReactionSystem)
-    D = Differential(rs.iv)
     eqs = essemble_drift(rs)
     noiseeqs = essemble_diffusion(rs)
-    SDESystem(eqs,noiseeqs,rs.iv,rs.dvs,rs.ps)
+    SDESystem(eqs,noiseeqs,rs.iv,rs.states,rs.ps,
+              name=rs.name,systems=convert.(SDESystem,rs.systems))
 end

test/reactions.jl test/reactionsystem.jltest/reactions.jl renamed to test/reactionsystem.jl

@@ -10,6 +10,7 @@ using SafeTestsets, Test
 @safetestset "SDESystem Test" begin include("sdesystem.jl") end
 @safetestset "NonlinearSystem Test" begin include("nonlinearsystem.jl") end
 @safetestset "OptimizationSystem Test" begin include("optimizationsystem.jl") end
+@safetestset "ReactionSystem Test" begin include("reactionsystem.jl") end
 @safetestset "Build Targets Test" begin include("build_targets.jl") end
 @safetestset "Domain Test" begin include("domains.jl") end
 @safetestset "Constraints Test" begin include("constraints.jl") end