add of @inbounds macro for perf

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")
 

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))

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
 
 """

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")

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"]