Merge pull request #49 from isaacsas/massaction-jumpproblems

Added conversions of jumps to MassActionJumps in JumpProblem

Author: ChrisRackauckas

REQUIRE

@@ -1,8 +1,8 @@
 julia 0.6
 DiffEqJump 4.2.0
-DiffEqBase 3.4.0
+DiffEqBase 3.11.0
 Compat 0.17.0
-DataStructures 0.4.6
+DataStructures 0.8.1
 Reexport 0.1.0
-SymEngine 0.3.0
+SymEngine 0.4.0
 MacroTools 0.4.0

src/DiffEqBiological.jl

@@ -14,6 +14,7 @@ using Compat
 
 include("reaction_network.jl")
 include("maketype.jl")
+include("massaction_jump_utils.jl")
 include("problem.jl")
 
 export @reaction_network, @reaction_func

src/massaction_jump_utils.jl

@@ -0,0 +1,79 @@
+# collection of utility functions for getting mass action jumps out of 
+# DiffEqBiological reaction network structures
+
+# return a map from species symbol to species index
+function species_to_indices(network)
+    specs = network.syms
+    Dict( specs[i] => i for i in eachindex(specs) )    
+end
+
+# return a map from reaction param symbol to rate index
+function rate_to_indices(network)
+    rates = network.params
+    Dict( rates[i] => i for i in eachindex(rates) )
+end
+
+# given a ReactionStruct and a species map construct a MassActionJump
+function make_majump(rs, specmap, ratemap, params)
+    reactant_stoich = Vector{Pair{Int,Int}}()
+    nsdict = Dict{Int,Int}()
+    for substrate in rs.substrates
+        specpair = specmap[substrate.reactant] => substrate.stoichiometry
+        push!(reactant_stoich, specpair)
+        specpair = specmap[substrate.reactant] => -substrate.stoichiometry
+        push!(nsdict, specpair)
+    end
+    sort!(reactant_stoich)
+
+    for product in rs.products
+        prodidx = specmap[product.reactant]
+        if haskey(nsdict, prodidx)
+            nsdict[prodidx] += product.stoichiometry
+        else
+            push!(nsdict, prodidx => product.stoichiometry)
+        end
+    end
+
+    net_stoich = Vector{Pair{Int,Int}}()
+    for stoich_map in sort(collect(nsdict))
+        push!(net_stoich, stoich_map)
+    end
+
+    if isempty(net_stoich)
+        error("Empty net stoichiometry vectors for mass action reactions are not allowed.")
+    end
+
+    if typeof(rs.rate_org) == Symbol
+        rateconst = params[ratemap[rs.rate_org]]
+    elseif typeof(rs.rate_org) == Expr
+        rateconst = eval(rs.rate_org)
+    elseif typeof(rs.rate_org) <: Number
+        rateconst = rs.rate_org
+    else
+        error("reaction_network.reactions.rate_org must have a type of Symbol, Expr, or Number.")
+    end
+
+    MassActionJump(Float64(rateconst), reactant_stoich, net_stoich)
+end
+
+# given a reaction network and species map, split the ConstantRateJumps and MassActionJumps
+function network_to_jumpset(rn, specmap, ratemap, params)
+
+    empty_majump = MassActionJump(0.0, [0=>1], [1=>1])
+    majumpvec    = Vector{typeof(empty_majump)}()
+    cjumpvec     = Vector{ConstantRateJump}()
+
+    for (i,rs) in enumerate(rn.reactions)        
+        if rs.is_pure_mass_action
+            push!(majumpvec, make_majump(rs, specmap, ratemap, params))
+        else
+            push!(cjumpvec, rn.jumps[i])
+        end
+    end
+
+    if isempty(majumpvec)
+        return JumpSet((), cjumpvec, nothing, nothing)
+    else
+        return JumpSet((), cjumpvec, nothing, majumpvec)
+    end
+end

src/problem.jl

@@ -3,11 +3,21 @@ DiffEqBase.SDEProblem(rn::AbstractReactionNetwork, u0::Union{AbstractArray, Numb
     SDEProblem(rn, rn.g, u0, args...;noise_rate_prototype=rn.p_matrix, kwargs...)
 
 ### JumpProblem ###
-function DiffEqJump.JumpProblem(prob,aggregator::Direct,rn::AbstractReactionNetwork; kwargs...)
+function DiffEqJump.JumpProblem(prob,aggregator,rn::AbstractReactionNetwork; kwargs...)
     if typeof(prob)<:DiscreteProblem && any(issubtype.(typeof.(rn.jumps),VariableRateJump))
         error("When using time dependant reaction rates a DiscreteProblem should not be used (try an ODEProblem). Also, use a continious solver.")
     end
-    JumpProblem(prob,aggregator::Direct,rn.jumps...;kwargs...)
+
+    # map from species symbol to index of species
+    spec_to_idx = species_to_indices(rn)
+
+    # map from parameter symbol to index of parameter in prob.p
+    param_to_idx = rate_to_indices(rn)
+
+    # get a JumpSet of the possible jumps
+    jset = network_to_jumpset(rn, spec_to_idx, param_to_idx, prob.p)
+
+    JumpProblem(prob, aggregator, jset; kwargs...)
 end
 
 ### SteadyStateProblem ###

test/mass_act_jump_tests.jl

@@ -0,0 +1,92 @@
+using DiffEqBiological, DiffEqJump, DiffEqBase, Base.Test
+
+dotestmean   = true
+doprintmeans = false
+reltol       = .01          # required test accuracy
+
+# run the given number of SSAs and return the mean
+function runSSAs(jump_prob, Nsims, idx)
+    Psamp = zeros(Int, Nsims)
+    for i in 1:Nsims
+        sol = solve(jump_prob, SSAStepper())
+        Psamp[i] = sol[idx,end]
+    end
+    mean(Psamp)
+end
+
+function execute_test(u0, tf, rates, rs, Nsims, expected_avg, idx, test_name)
+    prob = DiscreteProblem(u0, (0.0, tf), rates)
+    jump_prob = JumpProblem(prob, Direct(), rs)
+    avg_val = runSSAs(jump_prob, Nsims, idx)
+    
+    if dotestmean
+        if doprintmeans
+            println(test_name, ": mean = ", avg_val, ", act_mean = ", expected_avg)
+        end
+        @test abs(avg_val - expected_avg) < reltol * expected_avg
+    end
+    
+end
+
+# nonlinear reaction test
+Nsims = 32000
+tf = .01
+u0 = [200, 100, 150]
+expected_avg = 84.876015624999994
+rs = @reaction_network dtype begin
+    k1, 2A --> B
+    k2, B --> 2A 
+    k3, A + B --> C
+    k4, C --> A + B
+    k5, 3C --> 3A
+end k1 k2 k3 k4 k5
+rates = [1., 2., .5, .75, .25]
+execute_test(u0, tf, rates, rs, Nsims, expected_avg, 1, "Nonlinear rx test")
+
+
+# DNA repression model 
+rs = @reaction_network ptype begin
+    k1, DNA --> mRNA + DNA
+    k2, mRNA --> mRNA + P
+    k3, mRNA --> 0
+    k4, P --> 0
+    k5, DNA + P --> DNAR
+    k6, DNAR --> DNA + P
+end k1 k2 k3 k4 k5 k6
+Nsims        = 8000
+tf           = 1000.0
+u0           = [1,0,0,0]
+expected_avg = 5.926553750000000e+02
+rates = [.5, (20*log(2.)/120.), (log(2.)/120.), (log(2.)/600.), .025, 1.]
+execute_test(u0, tf, rates, rs, Nsims, expected_avg, 3, "DNA test")
+
+
+# simple constant production with degratation
+rs = @reaction_network pdtype begin
+    1000., 0 --> A
+    10, A --> 0
+end
+rates = nothing
+Nsims        = 16000
+tf           = 1.0
+u0           = [0]
+expected_avg = 1000./10*(1. - exp(-10*tf))
+execute_test(u0, tf, rates, rs, Nsims, expected_avg, 1, "Zero order test")
+
+
+# this is just a test to make sure a mixture of jumps actually runs without crashes
+network = @reaction_network rnType  begin
+    0.01, (X,Y,Z) --> 0
+    hill(X,3.,100.,-4), 0 --> Y
+    hill(Y,3.,100.,-4), 0 --> Z
+    hill(Z,4.5,100.,-4), 0 --> X
+    hill(X,2.,100.,6), 0 --> R
+    hill(Y,15.,100.,4)*0.002, R --> 0
+    20, 0 --> S
+    R*0.005, S --> SP
+    0.01, SP + SP --> SP2
+    0.05, SP2 --> 0
+end;
+prob = DiscreteProblem([200.,60.,120.,100.,50.,50.,50.], (0.,4000.))
+jump_prob = JumpProblem(prob, Direct(), network)
+sol = solve(jump_prob,SSAStepper());

test/runtests.jl

@@ -8,4 +8,5 @@ tic()
 @time @testset "Higher Order" begin include("higher_order_reactions.jl") end
 @time @testset "Additional Functions" begin include("func_test.jl") end
 @time @testset "Steady state solver" begin include("steady_state.jl") end
+@time @testset "Mass Action Jumps" begin include("mass_act_jump_tests.jl") end
 toc()