Setup ReactionSystem on AbstractSystem interface

Author: ChrisRackauckas

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

@@ -9,7 +9,8 @@ using SafeTestsets, Test
 @safetestset "Mass Matrix Test" begin include("mass_matrix.jl") end
 @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 "OptimizationSystem Test" begin include("optimizationsystem.jl") en
+@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