Merge pull request #46 from Gaussia/regular_jumps

Regular jumps

Author: ChrisRackauckas

src/maketype.jl

@@ -5,6 +5,7 @@ function maketype(name,
                   g,
                   g_func,
                   jumps,
+                  regular_jumps,
                   jump_rate_expr,
                   jump_affect_expr,
                   p_matrix,
@@ -22,6 +23,7 @@ function maketype(name,
         g::Function
         g_func::Vector{Any}
         jumps::Tuple{AbstractJump,Vararg{AbstractJump}}
+        regular_jumps::RegularJump
         jump_rate_expr::Tuple{Any,Vararg{Any}}
         jump_affect_expr::Tuple{Vector{Expr},Vararg{Vector{Expr}}}
         p_matrix::Array{Float64,2}
@@ -37,6 +39,7 @@ function maketype(name,
                   $(Expr(:kw,:g,g)),
                   $(Expr(:kw,:g_func,g_func)),
                   $(Expr(:kw,:jumps,jumps)),
+                  $(Expr(:kw,:regular_jumps,regular_jumps)),
                   $(Expr(:kw,:jump_rate_expr,jump_rate_expr)),
                   $(Expr(:kw,:jump_affect_expr,jump_affect_expr)),
                   $(Expr(:kw,:p_matrix,p_matrix)),
@@ -52,6 +55,7 @@ function maketype(name,
                       g,
                       g_func,
                       jumps,
+                      regular_jumps,
                       jump_rate_expr,
                       jump_affect_expr,
                       p_matrix,

src/reaction_network.jl

@@ -92,7 +92,7 @@ function coordinate(name, ex::Expr, p, scale_noise)
     g = make_func(g_expr, reactants, parameters)
     p_matrix = zeros(length(reactants), length(reactions))
 
-    (jump_rate_expr, jump_affect_expr, jumps) = get_jumps(reactions, reactants,parameters)
+    (jump_rate_expr, jump_affect_expr, jumps, regular_jumps) = get_jumps(reactions, reactants, parameters)
 
     f_rhs = [element.args[2] for element in f_expr]
     #symjac = Expr(:quote, calculate_jac(f_rhs, syms))
@@ -105,7 +105,7 @@ function coordinate(name, ex::Expr, p, scale_noise)
     f_funcs = [element.args[2] for element in f_expr]
     g_funcs = [element.args[2] for element in g_expr]
 
-    typeex,constructorex = maketype(name, f, f_funcs, f_symfuncs, g, g_funcs, jumps, Meta.quot(jump_rate_expr), Meta.quot(jump_affect_expr), p_matrix, syms; params=params, reactions=reactions)
+    typeex,constructorex = maketype(name, f, f_funcs, f_symfuncs, g, g_funcs, jumps, regular_jumps, Meta.quot(jump_rate_expr), Meta.quot(jump_affect_expr), p_matrix, syms; params=params, reactions=reactions)
 
     push!(exprs,typeex)
     push!(exprs,constructorex)
@@ -313,12 +313,15 @@ function get_jumps(reactions::Vector{ReactionStruct}, reactants::OrderedDict{Sym
     rates = Vector{Any}(length(reactions))
     affects = Vector{Vector{Expr}}(length(reactions))
     jumps = Expr(:tuple)
+    reg_rates = Expr(:block)
+    reg_c = Expr(:block)
     idx = 0
     for reaction in deepcopy(reactions)
         rates[idx += 1] = recursive_clean!(reaction.rate_SSA)
         affects[idx] = Vector{Expr}(0)
-        foreach(prod -> push!(affects[idx],:(@inbounds integrator.u[$(reactants[prod.reactant])] += $(prod.stoichiometry))), reaction.products)
-        foreach(sub -> push!(affects[idx],:(@inbounds integrator.u[$(reactants[sub.reactant])] -= $(sub.stoichiometry))), reaction.substrates)
+        reactant_set = union(getfield.(reaction.products, :reactant),getfield.(reaction.substrates, :reactant))
+        foreach(r -> push!(affects[idx],:(@inbounds integrator.u[$(reactants[r])] += $(get_stoch_diff(reaction,r)))), reactant_set)
+        syntax_rate = recursive_replace!(deepcopy(rates[idx]), (reactants,:internal_var___u), (parameters, :internal_var___p))
         #if reaction.is_pure_mass_action
         #    ma_sub_stoch = :(reactant_stoich = [[]])
         #    ma_stoch_change = :(reactant_stoich = [[]])
@@ -327,11 +330,14 @@ function get_jumps(reactions::Vector{ReactionStruct}, reactants::OrderedDict{Sym
         #    push!(jumps.args,:(MassActionJump($(reaction.rate_org),$(ma_sub_stoch),$(ma_stoch_change))))
         #else
             recursive_contains(:t,rates[idx]) ? push!(jumps.args,Expr(:call,:VariableRateJump)) : push!(jumps.args,Expr(:call,:ConstantRateJump))
-            push!(jumps.args[idx].args, :((internal_var___u,internal_var___p,t) -> $(recursive_replace!(deepcopy(rates[idx]), (reactants,:internal_var___u), (parameters, :internal_var___p)))))
+            push!(jumps.args[idx].args, :((internal_var___u,internal_var___p,t) -> $syntax_rate))
             push!(jumps.args[idx].args, :(integrator -> $(expr_arr_to_block(deepcopy(affects[idx])))))
         #end
+        push!(reg_rates.args,:(internal_var___out[$idx]=$syntax_rate))
+        foreach(r -> push!(reg_c.args,:(internal_var___dc[$(reactants[r]),$idx]=$(get_stoch_diff(reaction,r)))), reactant_set)
     end
-    return (Tuple(rates),Tuple(affects),jumps)
+    reg_jumps = :(RegularJump((internal_var___out,internal_var___u,internal_var___p,t)->$reg_rates,(internal_var___dc,internal_var___u,internal_var___p,t,internal_var___mark)->$reg_c,zeros($(length(reactants)),$(length(reactions)));constant_c=true))
+    return (Tuple(rates),Tuple(affects),jumps,reg_jumps)
 end
 
 #Recursively traverses an expression and removes things like X^1, 1*X. Will not actually have any affect on the expression when used as a function, but will make it much easier to look at it for debugging, as well as if it is transformed to LaTeX code.