From a4427ce9226f2d25a8e741160cb9cf691792b22d Mon Sep 17 00:00:00 2001
From: Samuel Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 8 May 2018 17:44:18 -0400
Subject: [PATCH] add diffusion example

---
 src/DiffEqProblemLibrary.jl  |  4 +++-
 src/jump_premade_problems.jl | 24 ++++++++++++++++++++++++
 2 files changed, 27 insertions(+), 1 deletion(-)

diff --git a/src/DiffEqProblemLibrary.jl b/src/DiffEqProblemLibrary.jl
index e01edfd..f7b45a4 100644
--- a/src/DiffEqProblemLibrary.jl
+++ b/src/DiffEqProblemLibrary.jl
@@ -51,6 +51,8 @@ export prob_jump_dnarepressor, prob_jump_constproduct, prob_jump_nonlinrxs,
 # examples mixing mass action and constant rate jumps
        prob_jump_osc_mixed_jumptypes,
 # examples used in published benchmarks / comparisions
-       prob_jump_multistate, prob_jump_twentygenes, prob_jump_dnadimer_repressor
+       prob_jump_multistate, prob_jump_twentygenes, prob_jump_dnadimer_repressor,
+# examples approximating diffusion by continuous time random walks
+       prob_jump_diffnetwork
 
 end # module
diff --git a/src/jump_premade_problems.jl b/src/jump_premade_problems.jl
index 4a135ff..df458b1 100644
--- a/src/jump_premade_problems.jl
+++ b/src/jump_premade_problems.jl
@@ -200,3 +200,27 @@ prob = DiscreteProblem(u0, (0.0, tf), rnpar)
 prob_jump_dnadimer_repressor = JumpProblemNetwork(rn, rnpar, tf, u0, prob, 
                                                 Dict("specs_names" => varlabels))
 
+
+"""
+    Continuous time random walk (i.e. diffusion approximation) example.
+    Here the network in the JumpProblemNetwork is a function that returns a 
+    network given the number of lattice sites.
+    u0 is a similar function that returns the initial condition vector.
+"""
+# diffusion model
+function getDiffNetwork(N)
+    diffnetwork = "@reaction_network dpldifftype begin\n"
+    for i in 1:(N-1)
+        diffnetwork *= "\t K, X$(i) --> X$(i+1)\n"
+        diffnetwork *= "\t K, X$(i+1) --> X$(i)\n"
+    end
+    diffnetwork *= "end K"
+    rs = eval( parse(diffnetwork) )
+    rs
+end
+params = (1.,)
+function getDiffu0(N)
+    10*ones(Int64, N)
+end
+tf = 10.
+prob_jump_diffnetwork = JumpProblemNetwork(getDiffNetwork, params, tf, getDiffu0, nothing, nothing)
\ No newline at end of file