Merge pull request #90 from isaacsas/use-discretefunctions

add DiscreteProblem overload for reaction_networks

Author: isaacsas

REQUIRE

@@ -1,6 +1,6 @@
 julia 1.0
 DiffEqJump 5.6.0
-DiffEqBase 3.11.0
+DiffEqBase 5.1.0
 Compat 0.17.0
 DataStructures 0.8.1
 Reexport 0.1.0

src/DiffEqBiological.jl

@@ -20,6 +20,6 @@ include("problem.jl")
 
 export @reaction_network, @reaction_func, @min_reaction_network
 export addodes!, addsdes!, addjumps!
-export ODEProblem, SDEProblem, JumpProblem, SteadyStateProblem
+export ODEProblem, SDEProblem, DiscreteProblem, JumpProblem, SteadyStateProblem
 
 end # module

src/problem.jl

@@ -10,6 +10,13 @@ function DiffEqBase.SDEProblem(rn::DiffEqBase.AbstractReactionNetwork, u0::Union
     SDEProblem(rn.sdefun, rn.g::Function, u0, args...;noise_rate_prototype=rn.p_matrix, kwargs...)
 end 
 
+### DiscreteProblem, passing through syms
+function DiffEqBase.DiscreteProblem(rn::DiffEqBase.AbstractReactionNetwork, u0, tspan::Tuple, p=nothing; kwargs...)
+    f = DiffEqBase.DISCRETE_INPLACE_DEFAULT
+    df = DiscreteFunction{true,typeof(f),Nothing,typeof(rn.syms)}(f,nothing,rn.syms)
+    DiscreteProblem(df, u0, tspan, p; kwargs...)
+end
+
 ### JumpProblem ###
 function build_jump_problem(prob, aggregator, rn, jumps, kwargs...)
     if typeof(prob)<:DiscreteProblem && any(x->typeof(x) <: VariableRateJump, jumps)

test/discreteproblem_test.jl

@@ -0,0 +1,66 @@
+using DiffEqBase, DiffEqBiological, DiffEqJump, Statistics
+
+dotestmean   = true
+doprintmeans = false
+reltol       = .01          # required test accuracy
+algs      = (Direct(),)
+
+# 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(prob, u0, tf, rates, rs, Nsims, expected_avg, idx, test_name)
+    for method in algs
+        jump_prob = JumpProblem(prob, method, rs, save_positions=(false,false))
+        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
+    end
+
+end
+
+# full macro
+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]
+prob = DiscreteProblem(u0, (0.0, tf), rates)
+execute_test(prob, u0, tf, rates, rs, Nsims, expected_avg, 1, "Nonlinear rx test (no syms)")
+
+prob = DiscreteProblem(rs, u0, (0.0, tf), rates)
+execute_test(prob, u0, tf, rates, rs, Nsims, expected_avg, 1, "Nonlinear rx test (syms)")
+
+# min_network macro
+rs = @min_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
+addjumps!(rs)
+prob = DiscreteProblem(u0, (0.0, tf), rates)
+execute_test(prob, u0, tf, rates, rs, Nsims, expected_avg, 1, "Nonlinear rx test, min network (no syms)")
+
+prob = DiscreteProblem(rs, u0, (0.0, tf), rates)
+execute_test(prob, u0, tf, rates, rs, Nsims, expected_avg, 1, "Nonlinear rx test, min network (syms)")
+

test/runtests.jl

@@ -1,6 +1,7 @@
 using DiffEqBiological
 using Test
 
+# full macro tests
 @time begin
   @time @testset "Model Macro" begin include("make_model_test.jl") end
   @time @testset "Gillespie Tests" begin include("gillespie.jl") end
@@ -11,6 +12,7 @@ using Test
   @time @testset "Mass Action Jumps" begin include("mass_act_jump_tests.jl") end
 end
 
+# min macro tests
 @time begin
   @time @testset "Model Macro (Min)" begin include("make_model_test_min.jl") end
   @time @testset "Gillespie Tests (Min)" begin include("gillespie_min.jl") end
@@ -19,4 +21,9 @@ end
   @time @testset "Additional Functions (Min)" begin include("func_test_min.jl") end
   @time @testset "Steady state solver (Min)" begin include("steady_state_min.jl") end
   @time @testset "Mass Action Jumps (Min)" begin include("mass_act_jump_tests_min.jl") end
+end
+
+# tests that handle both macros
+@time begin
+  @time @testset "Discrete Problem" begin include("discreteproblem_test.jl") end
 end