fix up SDESystem conversion

Author: ChrisRackauckas

src/ModelingToolkit.jl

@@ -91,6 +91,7 @@ include("systems/diffeqs/first_order_transform.jl")
 include("systems/diffeqs/modelingtoolkitize.jl")
 include("systems/nonlinear/nonlinear_system.jl")
 include("systems/pde/pdesystem.jl")
+include("systems/reaction/reactionsystem.jl")
 include("latexify_recipes.jl")
 include("build_function.jl")
 
@@ -99,6 +100,7 @@ export SDESystem, SDEFunction
 export NonlinearSystem
 export ode_order_lowering
 export PDESystem
+export Reaction, ReactionSystem
 export Differential, expand_derivatives, @derivatives
 export IntervalDomain, ProductDomain, ⊗, CircleDomain
 export Equation, ConstrainedEquation

src/systems/reaction/reactionsystem.jl

@@ -11,21 +11,53 @@ struct ReactionSystem
     ps
 end
 
-function Base.convert(::Type{<:ODESystem},rs::ReactionSystem)
+# 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]
 
     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)
+            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)
+            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]
+    @show size(eqs)
+
+    for (j,rx) in enumerate(rs.rxs)
+        for reactant in rx.reactants
+            i = findfirst(x->x.op == reactant.op,rs.dvs)
+            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)
-            eqs[i] = Equation(eqs[i].lhs,+(eqs[i].rhs,rx.reactants...))
+            eqs[i,j] += sqrt(rx.rate) * reaction_expr(rx.reactants)
         end
     end
+    eqs
+end
 
+function Base.convert(::Type{<:ODESystem},rs::ReactionSystem)
+    eqs = essemble_drift(rs)
     ODESystem(eqs,rs.iv,rs.dvs,rs.ps)
 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)
+end

test/reactions.jl

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