From 99ece473efd5e44b1f16f0628be3ec7f1e860dcf Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Wed, 24 Feb 2021 16:08:35 +0100
Subject: [PATCH] Big improvement of network models performance with SVector
---
bench/pkg/abstract_array.jl | 30 +++++++++----------
bench/pkg/euclidean_distance_repressilator.jl | 3 +-
bench/pkg/sim_repressilator.jl | 18 +++++++++++
core/MarkovProcesses.jl | 3 +-
core/network_model.jl | 6 ++--
models/ER.jl | 12 ++++----
models/SIR.jl | 10 +++----
models/SIR_tauleap.jl | 6 ++--
8 files changed, 53 insertions(+), 35 deletions(-)
create mode 100644 bench/pkg/sim_repressilator.jl
diff --git a/bench/pkg/abstract_array.jl b/bench/pkg/abstract_array.jl
index 2e163ee..13b0879 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 fd55221..2dacfd7 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 0000000..56fa166
--- /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 e848313..257be50 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 455d243..4a9468e 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 bd5d23c..193a4df 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 240ab49..7bc5336 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 01af738..e6743b5 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
--
GitLab