Merge pull request #346 from TorkelE/reaction_system_2_jump_system

Enable creation of JumpSystems from ReactionSystems

Author: ChrisRackauckas

src/systems/reaction/reactionsystem.jl

@@ -121,6 +121,49 @@ function assemble_diffusion(rs)
     eqs
 end
 
+function assemble_jumps(rs)
+    eqs = Vector{Union{ConstantRateJump, MassActionJump, VariableRateJump}}()
+
+    for rx in rs.eqs
+        rl = jumpratelaw(rx)
+        affect = Vector{Equation}()
+        for (spec,stoich) in rx.netstoich
+            push!(affect,var2op(spec) ~ var2op(spec) + stoich)
+        end
+        if any(isequal.(var2op(rs.iv),get_variables(rx.rate)))
+            push!(eqs,VariableRateJump(rl,affect))
+        elseif rx.only_use_rate || any([isequal(state,r_op) for state in rs.states, r_op in getfield.(get_variables(rx.rate),:op)])
+            push!(eqs,ConstantRateJump(rl,affect))
+        else
+            reactant_stoch = isempty(rx.substoich) ? [0 => 1] : Pair.(var2op.(getfield.(rx.substrates,:op)),rx.substoich)
+            net_stoch = map(p -> Pair(var2op(p[1]),p[2]),rx.netstoich)
+            push!(eqs,MassActionJump(rx.rate, reactant_stoch, net_stoch))
+        end
+    end
+    eqs
+end
+
+# Calculate the Jump rate law (like ODE, but uses X instead of X(t).
+# The former generates a "MethodError: objects of type Int64 are not callable" when trying to solve the problem.
+function jumpratelaw(rx)
+    @unpack rate, substrates, substoich, only_use_rate = rx
+    rl = deepcopy(rate)
+    foreach(op -> substitute_expr!(rl,op=>var2op(op.op)), get_variables(rx.rate))
+    if !only_use_rate
+        coef = one(eltype(substoich))
+        for (i,stoich) in enumerate(substoich)
+            coef *= factorial(stoich)
+            rl   *= isone(stoich) ? var2op(substrates[i].op) : var2op(substrates[i].op)^stoich
+        end
+        (!isone(coef)) && (rl /= coef)
+    end
+    rl
+end
+
+function var2op(var)
+    Operation(var,Vector{Expression}())
+end
+
 function Base.convert(::Type{<:ODESystem},rs::ReactionSystem)
     eqs = assemble_drift(rs)
     ODESystem(eqs,rs.iv,rs.states,rs.ps,name=rs.name,
@@ -133,3 +176,9 @@ function Base.convert(::Type{<:SDESystem},rs::ReactionSystem)
     SDESystem(eqs,noiseeqs,rs.iv,rs.states,rs.ps,
               name=rs.name,systems=convert.(SDESystem,rs.systems))
 end
+
+function Base.convert(::Type{<:JumpSystem},rs::ReactionSystem)
+    eqs = assemble_jumps(rs)
+    JumpSystem(eqs,rs.iv,rs.states,rs.ps,name=rs.name,
+              systems=convert.(JumpSystem,rs.systems))
+end