Skip to content
Snippets Groups Projects
Commit bae6e938 authored by Bentriou Mahmoud's avatar Bentriou Mahmoud
Browse files

change of model names

parent 300e9ec8
No related branches found
No related tags found
No related merge requests found
...@@ -4,14 +4,17 @@ using MarkovProcesses ...@@ -4,14 +4,17 @@ using MarkovProcesses
using Catalyst using Catalyst
using DiffEqJump using DiffEqJump
path_latex = "./"
# Bench 1: A long simulation of ER
load_model("ER") load_model("ER")
set_param!(ER, :k1, 0.2) set_param!(ER, :k1, 0.2)
set_param!(ER, :k2, 40.0) set_param!(ER, :k2, 40.0)
ER.buffer_size = 100 ER.buffer_size = 100
ER.estim_min_states = 8000 ER.estim_min_states = 8000
b_pkg = @benchmark simulate(ER) bench1_pkg = @benchmark simulate(ER)
@show minimum(b_pkg), mean(b_pkg), maximum(b_pkg)
rs = @reaction_network begin rs = @reaction_network begin
c1, S + E --> SE c1, S + E --> SE
...@@ -21,14 +24,78 @@ end c1 c2 c3 ...@@ -21,14 +24,78 @@ end c1 c2 c3
p = (0.2,40.0,1.0) # [c1,c2,c3] p = (0.2,40.0,1.0) # [c1,c2,c3]
tspan = (0., 100.) tspan = (0., 100.)
u0 = [100., 100., 0., 0.] # [S,E,SE,P] u0 = [100., 100., 0., 0.] # [S,E,SE,P]
# solve JumpProblem # solve JumpProblem
dprob = DiscreteProblem(rs, u0, tspan, p) dprob = DiscreteProblem(rs, u0, tspan, p)
jprob = JumpProblem(rs, dprob, Direct()) jprob = JumpProblem(rs, dprob, Direct())
jsol = solve(jprob, SSAStepper()) jsol = solve(jprob, SSAStepper())
b_catalyst = @benchmark solve(jprob, SSAStepper()) bench1_catalyst = @benchmark solve(jprob, SSAStepper())
@show minimum(b_catalyst), mean(b_catalyst), maximum(b_catalyst)
str_latex_bench1 = "
\\begin{tabular}{|c|c|c|c|c|}
\\hline
Bench 1 & Mean time (ms) & Max. time (ms) &
Min. time (ms) & \\begin{tabular}[c]{@{}c@{}}Mean\\\\Memory (KB)\\end{tabular} \\\\
\\hline
Package & $(round(time(mean(bench1_pkg))/1E6, digits=2)) & $(round(time(maximum(bench1_pkg))/1E6, digits=2)) &
$(round(time(minimum(bench1_pkg))/1E6, digits=2)) & $(round(memory(mean(bench1_pkg))/1024, digits=2)) \\\\
\\hline
Catalyst.jl & $(round(time(mean(bench1_catalyst))/1E6, digits=2)) & $(round(time(maximum(bench1_catalyst))/1E6, digits=2)) &
$(round(time(minimum(bench1_catalyst))/1E6, digits=2)) & $(round(memory(mean(bench1_catalyst))/1024, digits=2)) \\\\
\\hline
\\end{tabular}"
# Bench 2: Creation of model + a long simulation of ER
bench2_pkg = @benchmark begin
ER = @network_model begin
R1: (E+S => ES, k1*E*S)
R2: (ES => E+S, k2*ES)
R3: (ES => E+P, k3*ES)
end "ER"
set_param!(ER, :k1, 0.2)
set_param!(ER, :k2, 40.0)
ER.buffer_size = 100
ER.estim_min_states = 8000
simulate(ER)
end
bench2_catalyst = @benchmark begin
rs = @reaction_network begin
c1, S + E --> SE
c2, SE --> S + E
c3, SE --> P + E
end c1 c2 c3
p = (0.2,40.0,1.0) # [c1,c2,c3]
tspan = (0., 100.)
u0 = [100., 100., 0., 0.] # [S,E,SE,P]
# solve JumpProblem
dprob = DiscreteProblem(rs, u0, tspan, p)
jprob = JumpProblem(rs, dprob, Direct())
jsol = solve(jprob, SSAStepper())
end
str_latex_bench2 = "
\\begin{tabular}{|c|c|c|c|c|}
\\hline
Bench 1 & Mean time (ms) & Max. time (ms) &
Min. time (ms) & \\begin{tabular}[c]{@{}c@{}}Mean\\\\Memory (KB)\\end{tabular} \\\\
\\hline
Package & $(round(time(mean(bench2_pkg))/1E6, digits=2)) & $(round(time(maximum(bench2_pkg))/1E6, digits=2)) &
$(round(time(minimum(bench2_pkg))/1E6, digits=2)) & $(round(memory(mean(bench2_pkg))/1024, digits=2)) \\\\
\\hline
Catalyst.jl & $(round(time(mean(bench2_catalyst))/1E6, digits=2)) & $(round(time(maximum(bench2_catalyst))/1E6, digits=2)) &
$(round(time(minimum(bench2_catalyst))/1E6, digits=2)) & $(round(memory(mean(bench2_catalyst))/1024, digits=2)) \\\\
\\hline
\\end{tabular}"
@show minimum(bench1_pkg), mean(bench1_pkg), maximum(bench1_pkg)
@show minimum(bench1_catalyst), mean(bench1_catalyst), maximum(bench1_catalyst)
#plot(jsol,lw=2,title="Gillespie: Michaelis-Menten Enzyme Kinetics") open(path_latex * "bench1.tex", "w+") do io
write(io, str_latex_bench1)
end;
open(path_latex * "bench2.tex", "w+") do io
write(io, str_latex_bench2)
end;
...@@ -52,7 +52,7 @@ isabsorbing_ER(p::Vector{Float64},xn::Vector{Int}) = ...@@ -52,7 +52,7 @@ isabsorbing_ER(p::Vector{Float64},xn::Vector{Int}) =
@inbounds(p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3] === 0.0) @inbounds(p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3] === 0.0)
g_ER = [:P] g_ER = [:P]
ER = ContinuousTimeModel(d,k,dict_var_ER,dict_p_ER,l_tr_ER,p_ER,x0_ER,t0_ER,ER_f!,isabsorbing_ER; g=g_ER,name="ER pkg") ER = ContinuousTimeModel(d,k,dict_var_ER,dict_p_ER,l_tr_ER,p_ER,x0_ER,t0_ER,ER_f!,isabsorbing_ER; g=g_ER,name="ER SSA pkg")
function create_ER(new_p::Vector{Float64}) function create_ER(new_p::Vector{Float64})
ER_new = deepcopy(ER) ER_new = deepcopy(ER)
......
...@@ -50,7 +50,7 @@ end ...@@ -50,7 +50,7 @@ end
isabsorbing_SIR(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0 isabsorbing_SIR(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
g_SIR = [:I] g_SIR = [:I]
SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_SIR,p_SIR,x0_SIR,t0_SIR,SIR_f!,isabsorbing_SIR; g=g_SIR, name="SIR pkg") SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_SIR,p_SIR,x0_SIR,t0_SIR,SIR_f!,isabsorbing_SIR; g=g_SIR, name="SIR SSA pkg")
function create_SIR(new_p::Vector{Float64}) function create_SIR(new_p::Vector{Float64})
SIR_new = deepcopy(SIR) SIR_new = deepcopy(SIR)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment