# Diffusion Model

##### Samuel Isaacson, Chris Rackauckas
using DiffEqBase, DiffEqBiological, DiffEqJump, Plots, Statistics, DataFrames
gr()
fmt = :png

:png


# Model and example solutions

Here we implement a 1D continuous time random walk approximation of diffusion for $N$ lattice sites on $\left[0,1\right]$, with reflecting boundary conditions at $x=0$ and $x=1$.

N = 256
h = 1 / N
rn = @empty_reaction_network
function getDiffNetwork!(rn,N)
for i = 1:N
end
for i = 1:N
(i < N) && addreaction!(rn, :β, (Symbol(:u,i)=>1,), (Symbol(:u,i+1)=>1,))
(i > 1) && addreaction!(rn, :β, (Symbol(:u,i)=>1,), (Symbol(:u,i-1)=>1,))
end
rn
end
getDiffNetwork!(rn,N)
rnpar = [1/(h*h)]
u0 = 10*ones(Int64, N)
tf = .01

0.01

methods = (Direct(),DirectFW(),SortingDirect(),NRM(),DirectCR(),RSSA())
shortlabels = [string(leg)[12:end-2] for leg in methods]
prob    = prob = DiscreteProblem(u0, (0.0, tf), rnpar)
ploth   = plot(reuse=false)
for (i,method) in enumerate(methods)
println("Benchmarking method: ", method)
jump_prob = JumpProblem(prob, method, rn, save_positions=(false,false))
sol = solve(jump_prob, SSAStepper(), saveat=tf/1000.)
plot!(ploth,sol.t,sol[Int(N//2),:],label=shortlabels[i], format=fmt)
end

Benchmarking method: DiffEqJump.Direct()
Benchmarking method: DiffEqJump.DirectFW()
Benchmarking method: DiffEqJump.SortingDirect()
Benchmarking method: DiffEqJump.NRM()
Benchmarking method: DiffEqJump.DirectCR()

plot!(ploth, title="Population at middle lattice site", xlabel="time",format=fmt)


# Benchmarking performance of the methods

function run_benchmark!(t, jump_prob, stepper)
sol = solve(jump_prob, stepper)
@inbounds for i in 1:length(t)
t[i] = @elapsed (sol = solve(jump_prob, stepper))
end
end

run_benchmark! (generic function with 1 method)

nsims = 50
benchmarks = Vector{Vector{Float64}}()
for method in methods
jump_prob = JumpProblem(prob, method, rn, save_positions=(false,false))
stepper = SSAStepper()
t = Vector{Float64}(undef,nsims)
run_benchmark!(t, jump_prob, stepper)
push!(benchmarks, t)
end

medtimes = Vector{Float64}(undef,length(methods))
stdtimes = Vector{Float64}(undef,length(methods))
avgtimes = Vector{Float64}(undef,length(methods))
for i in 1:length(methods)
medtimes[i] = median(benchmarks[i])
avgtimes[i] = mean(benchmarks[i])
stdtimes[i] = std(benchmarks[i])
end

df = DataFrame(names=shortlabels,medtimes=medtimes,relmedtimes=(medtimes/medtimes[1]),
avgtimes=avgtimes, std=stdtimes, cv=stdtimes./avgtimes)


6 rows × 6 columns

namesmedtimesrelmedtimesavgtimesstdcv
StringFloat64Float64Float64Float64Float64
1Direct3.718031.03.718410.005263790.0014156
2DirectFW4.722711.270224.730080.02162150.00457106
3SortingDirect0.9153950.2462040.9151130.004693930.00512935
4NRM0.4145830.1115060.4153340.003052060.00734845
5DirectCR0.3413370.09180560.3424790.005082130.0148392

# Plotting

sa = [string(round(mt,digits=4),"s") for mt in df.medtimes]
bar(df.names,df.relmedtimes,legend=:false, fmt=fmt)
scatter!(df.names, .05 .+ df.relmedtimes, markeralpha=0, series_annotations=sa, fmt=fmt)
ylabel!("median relative to Direct")
title!("256 Site 1D Diffusion CTRW")

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


## Appendix

These benchmarks are a part of the DiffEqBenchmarks.jl repository, found at: https://github.com/JuliaDiffEq/DiffEqBenchmarks.jl

To locally run this tutorial, do the following commands:

using DiffEqBenchmarks
DiffEqBenchmarks.weave_file("Jumps","Diffusion_CTRW.jmd")

Computer Information:

Julia Version 1.3.0
Commit 46ce4d7933 (2019-11-26 06:09 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)
CPU: Intel(R) Core(TM) i7-9700K CPU @ 3.60GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-6.0.1 (ORCJIT, skylake)


Package Information:

Status: /home/julialab/.julia/dev/DiffEqBenchmarks/Project.toml
[28f2ccd6-bb30-5033-b560-165f7b14dc2f] ApproxFun 0.11.8
[a134a8b2-14d6-55f6-9291-3336d3ab0209] BlackBoxOptim 0.5.0
[a93c6f00-e57d-5684-b7b6-d8193f3e46c0] DataFrames 0.20.0
[2b5f629d-d688-5b77-993f-72d75c75574e] DiffEqBase 6.7.0
[eb300fae-53e8-50a0-950c-e21f52c2b7e0] DiffEqBiological 4.1.0
[f3b72e0c-5b89-59e1-b016-84e28bfd966d] DiffEqDevTools 2.16.0
[c894b116-72e5-5b58-be3c-e6d8d4ac2b12] DiffEqJump 6.3.0
[1130ab10-4a5a-5621-a13d-e4788d82bd4c] DiffEqParamEstim 1.10.0
[a077e3f3-b75c-5d7f-a0c6-6bc4c8ec64a9] DiffEqProblemLibrary 4.6.0
[ef61062a-5684-51dc-bb67-a0fcdec5c97d] DiffEqUncertainty 1.3.0
[0c46a032-eb83-5123-abaf-570d42b7fbaa] DifferentialEquations 6.9.0
[7f56f5a3-f504-529b-bc02-0b1fe5e64312] LSODA 0.6.1
[76087f3c-5699-56af-9a33-bf431cd00edd] NLopt 0.5.1
[c030b06c-0b6d-57c2-b091-7029874bd033] ODE 2.6.0
[54ca160b-1b9f-5127-a996-1867f4bc2a2c] ODEInterface 0.4.6
[09606e27-ecf5-54fc-bb29-004bd9f985bf] ODEInterfaceDiffEq 3.5.0
[1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] OrdinaryDiffEq 5.26.4
[65888b18-ceab-5e60-b2b9-181511a3b968] ParameterizedFunctions 4.2.1
[91a5bcdd-55d7-5caf-9e0b-520d859cae80] Plots 0.28.3
[b4db0fb7-de2a-5028-82bf-5021f5cfa881] ReactionNetworkImporters 0.1.5
[f2c3362d-daeb-58d1-803e-2bc74f2840b4] RecursiveFactorization 0.1.0
[9a3f8284-a2c9-5f02-9a11-845980a1fd5c] Random