diff --git a/bench/pkg/abstract_array.jl b/bench/pkg/abstract_array.jl index 2e163eeacad214bee3c34ae01238ff512a5f04c8..13b0879334020d33bc1ff6c565a1243710c70abc 100644 --- a/bench/pkg/abstract_array.jl +++ b/bench/pkg/abstract_array.jl @@ -3,25 +3,25 @@ using Profile using StaticArrays using BenchmarkTools -println("Cost of tests:") +println("Cost of 20 boolean tests:") f_test(a::Int) = a == 2 @btime for i = 1:20 f_test(3) end -# First container type +# Multiple dispatch abstract type SuperType end -for i = 1:10 +for i = 1:20 include_string(Main, "struct Type$i <: SuperType v::Symbol end") include_string(Main, "f(item::Type$i, a::Int) = a == $(i)") end vec_types = SuperType[] for i = 1:20 - rand_number_type = rand(1:10) + rand_number_type = rand(1:20) @eval item = $(Symbol("Type$(rand_number_type)"))(:test) push!(vec_types, item) end -println("First super type") +println("Multiple dispatch:") function read(vec_types::Vector{SuperType}, a::Int) for item in vec_types f(item, a) @@ -29,19 +29,19 @@ function read(vec_types::Vector{SuperType}, a::Int) end @btime read(vec_types, 3) -# Second container type +# Parametric struct struct AnotherSuperType{T<:Function} v::Symbol f::T end vec_others_types = AnotherSuperType[] for i = 1:20 - rand_number_type = rand(1:10) + rand_number_type = rand(1:20) include_string(Main, "f_other_$(i)(a::Int) = a == $(i)") sym_func = Symbol("f_other_$i") @eval push!(vec_others_types, AnotherSuperType(:test, $(sym_func))) end -println("Second super type") +println("Parametrized struct:") function read(vec_others_types::Vector{AnotherSuperType}, a::Int) for item in vec_others_types item.f(a) @@ -49,12 +49,12 @@ function read(vec_others_types::Vector{AnotherSuperType}, a::Int) end @btime read(vec_others_types, 3) -# With vectors -println("Transitions first:") +# With two vectors +println("Vectors:") vec_tr = Union{Nothing,Symbol}[] vec_func = Function[] for i = 1:20 - rand_number_type = rand(1:10) + rand_number_type = rand(1:20) include_string(Main, "f_vec_$(i)(a::Int) = a == $(i)") sym_func = Symbol("f_vec_$i") push!(vec_tr, :test) @@ -62,17 +62,17 @@ for i = 1:20 end function read(vec_tr::Vector{Union{Nothing,Symbol}}, vec_func::Vector{Function}, a::Int) for i = eachindex(vec_tr) - my_func = vec_func[i] - my_func(a) + vec_func[i](a) end end @btime read(vec_tr, vec_func, 3) +exit() println("Transitions second:") vec_tr = Union{Nothing,Symbol}[] vec_func = Function[] for i = 1:20 - rand_number_type = rand(1:10) + rand_number_type = rand(1:20) name_func = Symbol("f_vec2_$(i)") str_func = "$(name_func)(a::Int) = a == $(rand_number_type)" include_string(Main, str_func) @@ -84,7 +84,7 @@ end # How to store and read efficiently abstract types ? d = Dict{Symbol,Dict{Symbol,Vector{Float64}}}() -for i = 1:10 +for i = 1:20 sym_a = Symbol("a$i") sym_b = Symbol("a$i") d[sym_a] = Dict{Symbol,Vector{Float64}}() diff --git a/bench/pkg/euclidean_distance_repressilator.jl b/bench/pkg/euclidean_distance_repressilator.jl index fd55221349e0a06342f04a6e69673e107b0524f7..2dacfd789e942ada41c1327289f393a985f428c8 100644 --- a/bench/pkg/euclidean_distance_repressilator.jl +++ b/bench/pkg/euclidean_distance_repressilator.jl @@ -17,8 +17,7 @@ y_obs = vectorize(simulate(repressilator), :P1, tml_obs) println("Vectorize:") b_vectorize = @benchmark (σ = simulate($(repressilator)); euclidean_distance(σ, :P1, tml_obs, y_obs)) -@btime (σ = simulate($(repressilator))) -@btime (euclidean_distance(σ, :P1, tml_obs, y_obs)) +@btime (σ = simulate($(repressilator)); euclidean_distance(σ, :P1, tml_obs, y_obs)) @show minimum(b_vectorize), mean(b_vectorize), maximum(b_vectorize) println("Automaton with 1 loc") diff --git a/bench/pkg/sim_repressilator.jl b/bench/pkg/sim_repressilator.jl new file mode 100644 index 0000000000000000000000000000000000000000..56fa16611cdb13c5ef7caa5777576e3ff0f3fcde --- /dev/null +++ b/bench/pkg/sim_repressilator.jl @@ -0,0 +1,18 @@ + +using Statistics +using BenchmarkTools +using MarkovProcesses +import LinearAlgebra: dot +import Distributions: Uniform + +load_automaton("euclidean_distance_automaton") +load_model("repressilator") +tb = 210.0 +tml_obs = 0:10.0:210.0 +set_param!(repressilator, [:α, :β, :n, :α0], [200.0, 2.0, 2.0, 0.0]) +set_time_bound!(repressilator, tb) + +b_sim = @benchmark σ = simulate($(repressilator)) +@btime σ = simulate($(repressilator)) +@show mean(b_sim).time, mean(b_sim).memory + diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl index e848313a7be9671b38ee0dffe0eb9f9b9200dbf0..257be50be470f153d944f878680fb0467a245490 100644 --- a/core/MarkovProcesses.jl +++ b/core/MarkovProcesses.jl @@ -10,13 +10,14 @@ import Distributed: @everywhere, @distributed import Distributions: Product, Uniform, Normal import Distributions: Distribution, Univariate, Continuous, UnivariateDistribution, MultivariateDistribution, product_distribution -import StaticArrays: SVector +import StaticArrays: SVector, @SVector ## Exports export Distribution, Product, Uniform, Normal # Common types and constructors export Observations, AbstractTrajectory, Trajectory, SynchronizedTrajectory +export SVector, @SVector export Model, ContinuousTimeModel, SynchronizedModel, ParametricModel export VariableModel, ParameterModel, Transition, TransitionSet export LHA, StateLHA, Edge, Location, VariableAutomaton diff --git a/core/network_model.jl b/core/network_model.jl index 455d2436c24bb57482d96860c6a70356d71e1060..4a9468eec1d6fe086ac0f00f1a2b0aeb9ad56086 100644 --- a/core/network_model.jl +++ b/core/network_model.jl @@ -154,7 +154,7 @@ macro network_model(expr_network,expr_name...) symbol_func_f! = Symbol("$(basename_func)_f!") str_expr_model_f! = "function $(symbol_func_f!)(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition}, xn::Vector{Int}, tn::Float64, p::Vector{Float64})\n\t" # Computation of nu and propensity functions in f! - str_l_a = "l_a = (" + str_l_a = "l_a = SVector(" str_test_isabsorbing = "@inbounds(" l_nu = [zeros(Int, dim_state) for i = 1:nbr_reactions] for (i, expr_reaction) in enumerate(list_expr_reactions) @@ -205,8 +205,8 @@ macro network_model(expr_network,expr_name...) str_expr_model_f! *= "return nothing\n\t" str_expr_model_f! *= "end\n\t" # Computation of array of transitions - str_expr_model_f! *= "l_nu = (" * reduce(*, ["nu_$i, " for i = 1:nbr_reactions])[1:(end-2)] * ")\n\t" - str_expr_model_f! *= "l_sym_R = $(Tuple(transitions))\n\t" + str_expr_model_f! *= "l_nu = SVector(" * reduce(*, ["nu_$i, " for i = 1:nbr_reactions])[1:(end-2)] * ")\n\t" + str_expr_model_f! *= "l_sym_R = SVector$(Tuple(transitions))\n\t" # Simulation of the reaction str_expr_model_f! *= "u1 = rand()\n\t" str_expr_model_f! *= "u2 = rand()\n\t" diff --git a/models/ER.jl b/models/ER.jl index bd5d23c98be04b794c661bea570c47b85a4ffbbb..193a4df89394a407e0ae56cdd5bc8dc97e0a787b 100644 --- a/models/ER.jl +++ b/models/ER.jl @@ -12,17 +12,17 @@ t0_ER = 0.0 @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) + l_a = SVector(a1, a2, a3) asum = sum(l_a) if asum == 0.0 copyto!(xnplus1, xn) return nothing end - nu_1 = (-1, -1, 1, 0) - nu_2 = (1, 1, -1, 0) - nu_3 = (1, 0, -1, 1) - l_nu = (nu_1, nu_2, nu_3) - l_str_R = (:R1, :R2, :R3) + nu_1 = SVector(-1, -1, 1, 0) + nu_2 = SVector(1, 1, -1, 0) + nu_3 = SVector(1, 0, -1, 1) + l_nu = SVector(nu_1, nu_2, nu_3) + l_str_R = SVector(:R1, :R2, :R3) u1 = rand() u2 = rand() diff --git a/models/SIR.jl b/models/SIR.jl index 240ab4907ff315f9a13a5a881b83e5ca75696d75..7bc53369a05243dcc496450569ab0d3ed9cac93f 100644 --- a/models/SIR.jl +++ b/models/SIR.jl @@ -11,16 +11,16 @@ t0_SIR = 0.0 xn::Vector{Int}, tn::Float64, p::Vector{Float64}) @inbounds a1 = p[1] * xn[1] * xn[2] @inbounds a2 = p[2] * xn[2] - l_a = (a1, a2) + l_a = SVector(a1, a2) asum = sum(l_a) if asum == 0.0 copyto!(xnplus1, xn) return nothing end - nu_1 = (-1, 1, 0) - nu_2 = (0, -1, 1) - l_nu = (nu_1, nu_2) - l_str_R = (:R1,:R2) + nu_1 = SVector(-1, 1, 0) + nu_2 = SVector(0, -1, 1) + l_nu = SVector(nu_1, nu_2) + l_str_R = SVector(:R1,:R2) u1 = rand() u2 = rand() diff --git a/models/SIR_tauleap.jl b/models/SIR_tauleap.jl index 01af738bfaaf54fac570d8593a4765e3b434ccb3..e6743b5d27f0935eef71f63af27e3d28f9424b0b 100644 --- a/models/SIR_tauleap.jl +++ b/models/SIR_tauleap.jl @@ -14,15 +14,15 @@ t0_SIR_tauleap = 0.0 @inbounds tau = p[3] @inbounds a1 = p[1] * xn[1] * xn[2] @inbounds a2 = p[2] * xn[2] - l_a = (a1, a2) + l_a = SVector(a1, a2) asum = sum(l_a) if asum == 0.0 copyto!(xnplus1, xn) return nothing end # column-major order - nu_1 = (-1, 1, 0) - nu_2 = (0, -1, 1) + nu_1 = SVector(-1, 1, 0) + nu_2 = SVector(0, -1, 1) nbr_R1 = (a1 > 0.0) ? rand(Poisson(a1*tau)) : 0.0 nbr_R2 = (a2 > 0.0) ? rand(Poisson(a2*tau)) : 0.0 for i = 1:3