Merge pull request #50 from isaacsas/dep-graph

WIP: generate dependency graph for reaction networks

Author: ChrisRackauckas

src/massaction_jump_utils.jl

@@ -13,17 +13,24 @@ function rate_to_indices(network)
     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)
+# get substrate stoichiometry for a reaction 
+function get_substrate_stoich(rs, specmap)
     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)
+    end
+    sort!(reactant_stoich)
+    reactant_stoich
+end
+
+# get the net stoichiometry for a reaction
+function get_net_stoich(rs, specmap)
+    nsdict = Dict{Int,Int}()
+    for substrate in rs.substrates
         specpair = specmap[substrate.reactant] => -substrate.stoichiometry
         push!(nsdict, specpair)
     end
-    sort!(reactant_stoich)
 
     for product in rs.products
         prodidx = specmap[product.reactant]
@@ -36,9 +43,18 @@ function make_majump(rs, specmap, ratemap, params)
 
     net_stoich = Vector{Pair{Int,Int}}()
     for stoich_map in sort(collect(nsdict))
-        push!(net_stoich, stoich_map)
+        if stoich_map[2] != zero(Int)
+            push!(net_stoich, stoich_map)
+        end
     end
 
+    net_stoich
+end
+
+# given a ReactionStruct and a species map construct a MassActionJump
+function make_majump(rs, specmap, ratemap, params)
+    reactant_stoich = get_substrate_stoich(rs, specmap)
+    net_stoich      = get_net_stoich(rs, specmap)
     if isempty(net_stoich)
         error("Empty net stoichiometry vectors for mass action reactions are not allowed.")
     end
@@ -76,4 +92,81 @@ function network_to_jumpset(rn, specmap, ratemap, params)
     else
         return JumpSet((), cjumpvec, nothing, majumpvec)
     end
+end
+
+
+######################### dependency graph utilities #########################
+
+# map from a reaction index to the corresponding jump index
+# assumes all mass action jumps are ordered before constant rate jumps
+function rxidxs_to_jidxs_map(rn, num_majumps)
+    majidx = 1
+    cjidx  = num_majumps + 1
+    numrxs = length(rn.reactions)
+    rxi_to_ji = zeros(Int, numrxs)
+    for i = 1:numrxs 
+        if rn.reactions[i].is_pure_mass_action
+            rxi_to_ji[i] = majidx
+            majidx      += 1
+        else
+            rxi_to_ji[i] = cjidx
+            cjidx       += 1
+        end
+    end
+    rxi_to_ji
+end
+
+# map from species to Set of jumps depending on that species
+function spec_to_dep_jumps_map(rn, specmap, rxidxs_to_jidxs)
+
+    numrxs  = length(rn.reactions)
+    numspec = length(specmap)
+
+    # map from a species to jumps that depend on it
+    spec_to_dep_jumps = [Set{Int}() for n = 1:numspec]
+    for rx in 1:numrxs                    
+        for specsym in rn.reactions[rx].dependants
+            push!(spec_to_dep_jumps[specmap[specsym]], rxidxs_to_jidxs[rx])
+        end
+    end
+
+    spec_to_dep_jumps
+end
+
+# given a reaction network and species map, construct a jump dependency graph
+function depgraph_from_network(rn, specmap, jset)
+
+    numrxs          = length(rn.reactions)
+    numspec         = length(specmap)
+    num_majumps     = get_num_majumps(jset.massaction_jump)
+
+    # map reaction indices to jump indices
+    rxidxs_to_jidxs = rxidxs_to_jidxs_map(rn, num_majumps)
+
+    # map from species to jumps that depend on it
+    spec_to_dep_jumps = spec_to_dep_jumps_map(rn, specmap, rxidxs_to_jidxs)
+
+    # create map from a jump to jumps depending on it
+    dep_sets = [SortedSet{Int}() for n = 1:numrxs]
+    for rx in 1:numrxs        
+        jidx = rxidxs_to_jidxs[rx]
+
+        # get the net reaction stoichiometry
+        net_stoich = get_net_stoich(rn.reactions[rx], specmap)
+        
+        # rx changes spec, hence rxs depending on spec depend on rx
+        for (spec,stoch) in net_stoich
+            for dependent_jump in spec_to_dep_jumps[spec]
+                push!(dep_sets[jidx], dependent_jump)
+            end
+        end
+    end
+
+    # convert to Vectors of Vectors
+    dep_graph = Vector{Vector{Int}}(numrxs)
+    for jidx in 1:numrxs
+        dep_graph[jidx] = [dep for dep in dep_sets[jidx]]
+    end
+
+    dep_graph
 end

src/problem.jl

@@ -17,7 +17,14 @@ function DiffEqJump.JumpProblem(prob,aggregator,rn::AbstractReactionNetwork; kwa
     # 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...)
+    # # construct dependency graph if needed by aggregator
+    if needs_depgraph(aggregator)
+        dep_graph = depgraph_from_network(rn, spec_to_idx, jset)
+    else
+        dep_graph = nothing
+    end
+
+    JumpProblem(prob, aggregator, jset; dep_graph=dep_graph, kwargs...)
 end
 
 ### SteadyStateProblem ###

test/mass_act_jump_tests.jl

@@ -3,6 +3,7 @@ using DiffEqBiological, DiffEqJump, DiffEqBase, Base.Test
 dotestmean   = true
 doprintmeans = false
 reltol       = .01          # required test accuracy
+methods      = (Direct(), SortingDirect())
 
 # run the given number of SSAs and return the mean
 function runSSAs(jump_prob, Nsims, idx)
@@ -16,14 +17,17 @@ 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)
+
+    for method in methods
+        jump_prob = JumpProblem(prob, method, rs)
+        avg_val = runSSAs(jump_prob, Nsims, idx)
+        
+        if dotestmean
+            if doprintmeans
+                println(test_name, ", method = ", typeof(method), ", mean = ", avg_val, ", act_mean = ", expected_avg)
+            end
+            @test abs(avg_val - expected_avg) < reltol * expected_avg
         end
-        @test abs(avg_val - expected_avg) < reltol * expected_avg
     end
 
 end
@@ -55,12 +59,29 @@ rs = @reaction_network ptype begin
 end k1 k2 k3 k4 k5 k6
 Nsims        = 8000
 tf           = 1000.0
-u0           = [1,0,0,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")
 
 
+# DNA repression model, mix of jump types
+rs = @reaction_network ptype begin
+    k1*DNA, 0 --> mRNA
+    k2*mRNA, 0 --> 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           = [0,0,0,0]
+u0[findfirst(rs.syms, :DNA)] = 1
+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, findfirst(rs.syms, :P), "DNA mixed jump test")
+
 # simple constant production with degratation
 rs = @reaction_network pdtype begin
     1000., 0 --> A
@@ -88,5 +109,7 @@ network = @reaction_network rnType  begin
     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());
+for method in methods
+    jump_prob = JumpProblem(prob, method, network)
+    sol = solve(jump_prob,SSAStepper());
+end