diff --git a/bench/pkg/catalyst.jl b/bench/pkg/catalyst.jl
index 1f81b00b92ac430ebed6cd5c49393520eb6c2d34..cac5bb6cec921e7df2e965f64e1d6d3f3848a455 100644
--- a/bench/pkg/catalyst.jl
+++ b/bench/pkg/catalyst.jl
@@ -11,6 +11,7 @@ ER.buffer_size = 100
ER.estim_min_states = 8000
b_pkg = @benchmark simulate(ER)
+@show minimum(b_pkg), mean(b_pkg), maximum(b_pkg)
rs = @reaction_network begin
c1, S + E --> SE
@@ -27,6 +28,7 @@ jprob = JumpProblem(rs, dprob, Direct())
jsol = solve(jprob, SSAStepper())
b_catalyst = @benchmark solve(jprob, SSAStepper())
+@show minimum(b_catalyst), mean(b_catalyst), maximum(b_catalyst)
#plot(jsol,lw=2,title="Gillespie: Michaelis-Menten Enzyme Kinetics")
diff --git a/core/biochemical_network.jl b/core/biochemical_network.jl
index 0866b412ee7e70530cfae8615de7030abfbdda6c..3ae5b5f32cbef35f877b7d1f0c017abe71e04d87 100644
--- a/core/biochemical_network.jl
+++ b/core/biochemical_network.jl
@@ -108,13 +108,13 @@ macro biochemical_network(expr_network,expr_name...)
expr_model_f! = "function $(basename_func)_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}}, xn::Vector{Int}, tn::Float64, p::Vector{Float64})\n\t"
# Computation of nu and propensity functions in f!
str_l_a = "l_a = ("
- str_test_isabsorbing = "("
+ str_test_isabsorbing = "@inbounds("
l_nu = [zeros(Int, dim_state) for i = 1:nbr_reactions]
for (i, expr_reaction) in enumerate(list_expr_reactions)
local isreaction = @capture(expr_reaction, TR_: (reactants_ => products_, propensity_))
# Writing of propensities function
str_propensity = get_str_propensity(propensity, dict_species, dict_params)
- expr_model_f! *= "a$(i) = " * str_propensity * "\n\t"
+ expr_model_f! *= "@inbounds a$(i) = " * str_propensity * "\n\t"
# Anticipating the write of the function isabsorbing
str_test_isabsorbing *= str_propensity * "+"
# Update the nu of the i-th reaction
@@ -164,17 +164,16 @@ macro biochemical_network(expr_network,expr_name...)
expr_model_f! *= "reaction = i\n\t\t\t"
expr_model_f! *= "break\n\t\t"
expr_model_f! *= "end\n\t\t"
- expr_model_f! *= "b_inf += l_a[i]\n\t\t"
- expr_model_f! *= "b_sup += l_a[i+1]\n\t"
+ expr_model_f! *= "@inbounds b_inf += l_a[i]\n\t\t"
+ expr_model_f! *= "@inbounds b_sup += l_a[i+1]\n\t"
expr_model_f! *= "end\n\t"
expr_model_f! *= "nu = l_nu[reaction]\n\t"
expr_model_f! *= "for i = 1:$(dim_state)\n\t\t"
- expr_model_f! *= "xnplus1[i] = xn[i]+nu[i]\n\t"
+ expr_model_f! *= "@inbounds xnplus1[i] = xn[i]+nu[i]\n\t"
expr_model_f! *= "end\n\t"
- expr_model_f! *= "l_t[1] = tn + tau\n\t"
- expr_model_f! *= "l_tr[1] = l_str_R[reaction]\n"
+ expr_model_f! *= "@inbounds l_t[1] = tn + tau\n\t"
+ expr_model_f! *= "@inbounds l_tr[1] = l_str_R[reaction]\n"
expr_model_f! *= "end\n"
-
expr_model_isabsorbing = "isabsorbing_$(basename_func)(p::Vector{Float64},xn::Vector{Int}) = $(str_test_isabsorbing) === 0.0"
model_f! = eval(Meta.parse(expr_model_f!))
model_isabsorbing = eval(Meta.parse(expr_model_isabsorbing))
diff --git a/core/model.jl b/core/model.jl
index 13809b739633b01e032150bf402dc3656142e653..be93b4d51f91555d1d5509d8fbad5b34dc223494 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -4,24 +4,23 @@ import Distributions: insupport, pdf
function _resize_trajectory!(values::Vector{Vector{Int}}, times::Vector{Float64},
transitions::Vector{Transition}, size::Int)
- for i = eachindex(values) resize!(values[i], size) end
+ for i = eachindex(values) resize!(@inbounds(values[i]), size) end
resize!(times, size)
resize!(transitions, size)
end
-
function _finish_bounded_trajectory!(values::Vector{Vector{Int}}, times::Vector{Float64},
transitions::Vector{Transition}, time_bound::Float64)
- for i = eachindex(values) push!(values[i], values[i][end]) end
+ for i = eachindex(values) push!(@inbounds(values[i]), @inbounds(values[i][end])) end
push!(times, time_bound)
push!(transitions, nothing)
end
function _update_values!(values::Vector{Vector{Int}}, times::Vector{Float64}, transitions::Vector{Transition},
xn::Vector{Int}, tn::Float64, tr_n::Transition, idx::Int)
- for k = eachindex(values) values[k][idx] = xn[k] end
- times[idx] = tn
- transitions[idx] = tr_n
+ for k = eachindex(values) @inbounds(values[k][idx] = xn[k]) end
+ @inbounds(times[idx] = tn)
+ @inbounds(transitions[idx] = tr_n)
end
"""
diff --git a/models/ER.jl b/models/ER.jl
index 2af5b5f409ec0155b6aa095a34d530e105335b4c..f41a80d25e0417f72143d43226d17bc9efed7ae2 100644
--- a/models/ER.jl
+++ b/models/ER.jl
@@ -11,9 +11,9 @@ x0_ER = [100, 100, 0, 0]
t0_ER = 0.0
function ER_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}},
xn::Vector{Int}, tn::Float64, p::Vector{Float64})
- a1 = p[1] * xn[1] * xn[2]
- a2 = p[2] * xn[3]
- a3 = p[3] * xn[3]
+ @inbounds a1 = p[1] * xn[1] * xn[2]
+ @inbounds a2 = p[2] * xn[3]
+ @inbounds a3 = p[3] * xn[3]
l_a = (a1, a2, a3)
asum = sum(l_a)
if asum == 0.0
@@ -37,19 +37,19 @@ function ER_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{No
reaction = i
break
end
- b_inf += l_a[i]
- b_sup += l_a[i+1]
+ @inbounds b_inf += l_a[i]
+ @inbounds b_sup += l_a[i+1]
end
nu = l_nu[reaction]
for i = 1:4
- xnplus1[i] = xn[i]+nu[i]
+ @inbounds xnplus1[i] = xn[i]+nu[i]
end
- l_t[1] = tn + tau
- l_tr[1] = l_str_R[reaction]
+ @inbounds l_t[1] = tn + tau
+ @inbounds l_tr[1] = l_str_R[reaction]
end
isabsorbing_ER(p::Vector{Float64},xn::Vector{Int}) =
- (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"]
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")
diff --git a/models/SIR.jl b/models/SIR.jl
index 592835f1d44553fa8b858ba686ee59fbbdfe6876..f53db51dae007fbd34476c6c1af73aeeef3717e5 100644
--- a/models/SIR.jl
+++ b/models/SIR.jl
@@ -11,8 +11,8 @@ x0_SIR = [95, 5, 0]
t0_SIR = 0.0
function SIR_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}},
xn::Vector{Int}, tn::Float64, p::Vector{Float64})
- a1 = p[1] * xn[1] * xn[2]
- a2 = p[2] * xn[2]
+ @inbounds a1 = p[1] * xn[1] * xn[2]
+ @inbounds a2 = p[2] * xn[2]
l_a = (a1, a2)
asum = sum(l_a)
if asum == 0.0
@@ -36,16 +36,16 @@ function SIR_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{N
reaction = i
break
end
- b_inf += l_a[i]
- b_sup += l_a[i+1]
+ @inbounds b_inf += l_a[i]
+ @inbounds b_sup += l_a[i+1]
end
nu = l_nu[reaction]
for i = 1:3
- xnplus1[i] = xn[i]+nu[i]
+ @inbounds xnplus1[i] = xn[i]+nu[i]
end
- l_t[1] = tn + tau
- l_tr[1] = l_str_R[reaction]
+ @inbounds l_t[1] = tn + tau
+ @inbounds l_tr[1] = l_str_R[reaction]
end
isabsorbing_SIR(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
g_SIR = ["I"]