Commit cefe75b7 authored by Bentriou Mahmoud's avatar Bentriou Mahmoud
Browse files

Change of different names as the Continuous Time model object.

Major fix of type unstabilities which improved performances.
For now on, we observe for pygmalion benchmarks:
+2x speed simulation compared to simulation of models.
100x speed in terms of reading values (no allocations).
parent 8e760732
......@@ -31,6 +31,13 @@ MarkovProcesses.load_model("ER")
ER.time_bound = 20.0
@timev MarkovProcesses.simulate(ER)
println("Default buffer size=10")
b1 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate($ER) end
b2 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end
@show minimum(b1), mean(b1), maximum(b1)
@show minimum(b2), mean(b2), maximum(b2)
println("Buffer size 100")
ER.buffer_size = 100
b1 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate($ER) end
b2 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end
@show minimum(b1), mean(b1), maximum(b1)
......
......@@ -16,7 +16,7 @@ u = Control(10.0)
tml = 1:400
g_all = create_observation_function([ObserverModel(str_oml, tml)])
so = pygmalion.simulate(f, g_all, x0, u, p_true; on = nothing, full_timeline = true)
function read_trajectory_v1(so::SystemObservation)
function read_trajectory_pyg_v1(so::SystemObservation)
vals = to_vec(so, "P")
n_states = length(vals)
res = 0.0
......@@ -25,7 +25,7 @@ function read_trajectory_v1(so::SystemObservation)
end
return res
end
function read_trajectory_v2(so::SystemObservation)
function read_trajectory_pyg_v2(so::SystemObservation)
idx_P = om_findfirst("P", so.oml)
n_states = length(so.otll[idx_P])
res = 0.0
......@@ -35,11 +35,11 @@ function read_trajectory_v2(so::SystemObservation)
return res
end
# Bench
@timev read_trajectory_v1(so)
b1_pyg = @benchmark read_trajectory_v1($so)
@timev read_trajectory_pyg_v1(so)
b1_pyg = @benchmark read_trajectory_pyg_v1($so)
@show minimum(b1_pyg), mean(b1_pyg), maximum(b1_pyg)
@timev read_trajectory_v2(so)
b2_pyg = @benchmark read_trajectory_v2($so)
@timev read_trajectory_pyg_v2(so)
b2_pyg = @benchmark read_trajectory_pyg_v2($so)
@show minimum(b2_pyg), mean(b2_pyg), maximum(b2_pyg)
println("MarkovProcesses:")
......@@ -47,16 +47,28 @@ println("MarkovProcesses:")
MarkovProcesses.load_model("ER")
ER.time_bound = 10.0
σ = MarkovProcesses.simulate(ER)
function read_trajectory(σ::AbstractTrajectory)
function read_trajectory_v1(σ::AbstractTrajectory)
n_states = length_states(σ)
res = 0
for i = 1:n_states
res += (σ["P"])[i]
res += σ["P",i]
end
return res
end
function read_trajectory_v2(σ::AbstractTrajectory)
n_states = length_states(σ)
vals = σ["P"]
res = 0
for i = 1:n_states
res += vals[i]
end
return res
end
# Bench
@timev read_trajectory(σ)
b1 = @benchmark read_trajectory($σ)
@timev read_trajectory_v1(σ)
b1 = @benchmark read_trajectory_v1($σ)
@show minimum(b1), mean(b1), maximum(b1)
@timev read_trajectory_v2(σ)
b2 = @benchmark read_trajectory_v2($σ)
@show minimum(b2), mean(b2), maximum(b2)
module MarkovProcesses
import Base: +, -, getfield, getindex
import Base: +, -, getfield, getindex, ImmutableDict
export Model, ContinuousTimeModel, DiscreteTimeModel
export simulate, set_param!, get_param, set_observed_var!
......@@ -10,7 +10,7 @@ include("model.jl")
export Observations, AbstractTrajectory, Trajectory
export +, -, δ, dist_lp
export get_obs_var, length_states, length_obs_var, is_bounded
export get_obs_var, length_states, length_obs_var, is_bounded, times, transitions
include("trajectory.jl")
end
......
......@@ -2,15 +2,14 @@
import StaticArrays: SVector
abstract type Model end
abstract type ContinuousTimeModel <: Model end
abstract type DiscreteTimeModel <: Model end
mutable struct CTMC <: ContinuousTimeModel
mutable struct ContinuousTimeModel <: Model
d::Int # state space dim
k::Int # parameter space dim
map_var_idx::Dict # maps str to full state space
_map_obs_var_idx::Dict # maps str to observed state space
map_param_idx::Dict # maps str in parameter space
map_var_idx::Dict{String,Int} # maps str to full state space
_map_obs_var_idx::Dict{String,Int} # maps str to observed state space
map_param_idx::Dict{String,Int} # maps str in parameter space
l_name_transitions::Vector{String}
p::Vector{Float64}
x0::Vector{Int}
......@@ -23,7 +22,7 @@ mutable struct CTMC <: ContinuousTimeModel
buffer_size::Int
end
function CTMC(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_name_transitions::Vector{String},
function ContinuousTimeModel(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_name_transitions::Vector{String},
p::Vector{Float64}, x0::Vector{Int}, t0::Float64,
f!::Function, is_absorbing::Function;
g::Vector{String} = keys(map_var_idx), time_bound::Float64 = Inf, buffer_size::Int = 10)
......@@ -42,7 +41,7 @@ function CTMC(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_name_tra
@warn "You have possibly redefined a function Model.is_absorbing used in a previously instantiated model."
end
return CTMC(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_name_transitions, p, x0, t0, f!, g, _g_idx, is_absorbing, time_bound, buffer_size)
return ContinuousTimeModel(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_name_transitions, p, x0, t0, f!, g, _g_idx, is_absorbing, time_bound, buffer_size)
end
function simulate(m::ContinuousTimeModel)
......@@ -53,27 +52,27 @@ function simulate(m::ContinuousTimeModel)
transitions = Union{String,Nothing}[nothing]
# values at time n
n = 0
xn = @view m.x0[:] # View for type stability
xn = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
tn = m.t0
# at time n+1
mat_x = zeros(Int, m.buffer_size, m.d)
l_t = zeros(Float64, m.buffer_size)
l_tr = Vector{String}(undef, m.buffer_size)
is_absorbing = m.is_absorbing(m.p,xn)::Bool
is_absorbing::Bool = m.is_absorbing(m.p,xn)
while !is_absorbing && (tn <= m.time_bound)
i = 0
while i < m.buffer_size && !is_absorbing && (tn <= m.time_bound)
i += 1
m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
xn = @view mat_x[i,:]
xn = view(mat_x, i, :)
tn = l_t[i]
is_absorbing = m.is_absorbing(m.p,xn)::Bool
is_absorbing = m.is_absorbing(m.p,xn)
end
full_values = vcat(full_values, @view mat_x[1:i,:])
append!(times, @view l_t[1:i])
append!(transitions, @view l_tr[1:i])
full_values = vcat(full_values, view(mat_x, 1:i, :))
append!(times, view(l_t, 1:i))
append!(transitions, view(l_tr, 1:i))
n += i
is_absorbing = m.is_absorbing(m.p,xn)::Bool
is_absorbing = m.is_absorbing(m.p,xn)
end
if is_bounded(m)
if times[end] > m.time_bound
......@@ -86,7 +85,7 @@ function simulate(m::ContinuousTimeModel)
push!(transitions, nothing)
end
end
values = @view full_values[:,m._g_idx]
values = view(full_values, :, m._g_idx)
return Trajectory(m, values, times, transitions)
end
......
......@@ -74,11 +74,11 @@ function dist_lp(σ1::AbstractTrajectory, σ2::AbstractTrajectory, var::String;
if !is_bounded(σ1) || !is_bounded(σ2)
@warn "Lp distance computed on unbounded trajectories. Result should be wrong"
end
return dist_lp(σ1[var], σ1["times"], σ2[var], σ2["times"]; verbose = false, p = p)
return dist_lp(σ1[var], times(σ1), σ2[var], times(σ2); verbose = false, p = p)
end
# Distance function. Vectorized version
function dist_lp(x_obs::SubArray{Int,1}, t_x::SubArray{Float64,1}, y_obs::SubArray{Int,1}, t_y::SubArray{Float64,1};
function dist_lp(x_obs::SubArray{Int,1}, t_x::Vector{Float64}, y_obs::SubArray{Int,1}, t_y::Vector{Float64};
verbose::Bool = false, p::Int = 1)
current_y_obs = y_obs[1]
current_t_y = t_y[2]
......@@ -179,32 +179,23 @@ is_bounded(σ::AbstractTrajectory) = σ.transitions[end] == nothing
# Access to trajectory values
get_var_values(σ::AbstractTrajectory, var::String) =
@view σ.values[:,(σ.m)._map_obs_var_idx[var]]
get_state(σ::AbstractTrajectory, idx::Int) = @view σ.values[idx,:]
view(σ.values, :, (σ.m)._map_obs_var_idx[var])
get_state(σ::AbstractTrajectory, idx::Int) =
view(σ.values, idx, :)
get_value(σ::AbstractTrajectory, var::String, idx::Int) =
σ.values[idx,(σ.m)._map_obs_var_idx[var]]
function δ(σ1::AbstractTrajectory,t::Float64) end
# Get var values ["I"]
function getindex(σ::AbstractTrajectory, var::String)
if var == "times"
return @view σ.times[:]
elseif var == "transitions"
return @view σ.transitions[:]
else
return get_var_values(σ, var)
end
end
σ.values[idx,(σ.m)._map_obs_var_idx[var]]
δ(σ::AbstractTrajectory,t::Float64) = true
states(σ::AbstractTrajectory) = σ.values
times(σ::AbstractTrajectory) = σ.times
transitions(σ::AbstractTrajectory) = σ.transitions
# Get i-th state [i]
getindex(σ::AbstractTrajectory, idx::Int) = get_state(σ, idx)
# Get i-th value of var ["I", idx]
function getindex(σ::AbstractTrajectory, var::String, idx::Int)
if var == "times"
return σ.times[idx]
elseif var == "transitions"
return σ.transitions[idx]
else
return get_value(σ, var, idx)
end
end
getindex(σ::AbstractTrajectory, var::String, idx::Int) = get_value(σ, var, idx)
# Get the path of a variable ["I"]
getindex(σ::AbstractTrajectory, var::String) = get_var_values(σ, var)
......@@ -45,6 +45,6 @@ is_absorbing_ER(p::Vector{Float64},xn::SubArray{Int,1}) =
(p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0
g = ["P"]
ER = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_f!,is_absorbing_ER; g=g)
ER = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_f!,is_absorbing_ER; g=g)
export ER
......@@ -34,7 +34,7 @@ function SIR_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String},
b_sup += l_a[i+1]
end
nu = @view l_nu[:,reaction] # macro for avoiding a copy
nu = view(l_nu, :, reaction) # macro for avoiding a copy
for i = 1:3
mat_x[idx,i] = xn[i]+nu[i]
end
......@@ -44,6 +44,6 @@ end
is_absorbing_SIR(p::Vector{Float64}, xn::SubArray{Int,1}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
g = ["I"]
SIR = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_f!,is_absorbing_SIR; g=g)
SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_f!,is_absorbing_SIR; g=g)
export SIR
......@@ -46,6 +46,6 @@ is_absorbing_ER_col(p::Vector{Float64},xn::AbstractVector{Int}) =
(p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0
g = ["P"]
ER_col = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_col_f!,is_absorbing_ER_col; g=g)
ER_col = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_col_f!,is_absorbing_ER_col; g=g)
export ER_col
......@@ -45,6 +45,6 @@ is_absorbing_ER_col_buffer(p::Vector{Float64},xn::AbstractVector{Int}) =
(p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0
g = ["P"]
ER_col_buffer = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_col_buffer_f!,is_absorbing_ER_col_buffer; g=g)
ER_col_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_col_buffer_f!,is_absorbing_ER_col_buffer; g=g)
export ER_col_buffer
......@@ -45,6 +45,6 @@ is_absorbing_ER_row_buffer(p::Vector{Float64},xn::AbstractVector{Int}) =
(p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0
g = ["P"]
ER_row_buffer = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_row_buffer_f!,is_absorbing_ER_row_buffer; g=g)
ER_row_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_row_buffer_f!,is_absorbing_ER_row_buffer; g=g)
export ER_row_buffer
......@@ -44,6 +44,6 @@ end
is_absorbing_SIR_col(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
g = ["I"]
SIR_col = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_f!,is_absorbing_SIR_col; g=g)
SIR_col = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_f!,is_absorbing_SIR_col; g=g)
export SIR_col
......@@ -44,6 +44,6 @@ end
is_absorbing_SIR_col_buffer(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
g = ["I"]
SIR_col_buffer = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_buffer_f!,is_absorbing_SIR_col_buffer; g=g)
SIR_col_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_buffer_f!,is_absorbing_SIR_col_buffer; g=g)
export SIR_col_buffer
......@@ -44,6 +44,6 @@ end
is_absorbing_SIR_row_buffer(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
g = ["I"]
SIR_row_buffer = CTMC(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_row_buffer_f!,is_absorbing_SIR_row_buffer; g=g)
SIR_row_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_row_buffer_f!,is_absorbing_SIR_row_buffer; g=g)
export SIR_row_buffer
......@@ -23,9 +23,9 @@ for p = 1:2
l_tr = Vector{Nothing}(nothing, length(x_obs))
σ1 = Trajectory(SIR, values, t_x, l_tr)
test_1 = isapprox(dist_lp(@view(x_obs[:]), @view(t_x[:]), @view(y_obs[:]), @view(t_y[:]); p=p), res; atol = 1E-10)
test_1 = isapprox(dist_lp(@view(x_obs[:]), t_x, @view(y_obs[:]), t_y; p=p), res; atol = 1E-10)
test_1 = test_1 && isapprox(dist_lp(σ1,σ2; p=p), res; atol = 1E-10)
test_1_bis = isapprox(dist_lp(@view(y_obs[:]), @view(t_y[:]), @view(x_obs[:]), @view(t_x[:]); p=p), res; atol = 1E-10)
test_1_bis = isapprox(dist_lp(@view(y_obs[:]), t_y, @view(x_obs[:]), t_x; p=p), res; atol = 1E-10)
test_1_bis = test_1_bis && isapprox(dist_lp(σ2,σ1; p=p), res; atol = 1E-10)
f_x(t::Real) = MarkovProcesses._f_step(x_obs, t_x, t)
......
......@@ -26,9 +26,9 @@ for p = 1:2
int_riemann = MarkovProcesses._riemann_sum(diff_f, 0.0, 20.0, 1E-5)
int_riemann = int_riemann^(1/p)
res1 = dist_lp(@view(x_obs[:]), @view(t_x[:]), @view(y_obs[:]), @view(t_y[:]); p=p)
res1 = dist_lp(@view(x_obs[:]), t_x, @view(y_obs[:]), t_y; p=p)
res2 = dist_lp(σ1,σ2; p=p)
res1_bis = dist_lp(@view(y_obs[:]), @view(t_y[:]), @view(x_obs[:]), @view(t_x[:]); p=p)
res1_bis = dist_lp(@view(y_obs[:]), t_y, @view(x_obs[:]), t_x; p=p)
res2_bis = dist_lp(σ2,σ1; p=p)
test_1 = isapprox(res1, int_riemann; atol = 1E-3)
test_1 = test_1 && isapprox(res2, int_riemann; atol = 1E-3)
......
......@@ -21,7 +21,7 @@ values[:,1] = y_obs
l_tr = Vector{Nothing}(nothing, length(y_obs))
σ2 = Trajectory(SIR, values, t_y, l_tr)
test_1 = dist_lp(@view(x_obs[:]), @view(t_x[:]), @view(y_obs[:]), @view(t_y[:]); p=1) == 6.4
test_1 = dist_lp(@view(x_obs[:]), t_x, @view(y_obs[:]), t_y; p=1) == 6.4
test_1 = test_1 && dist_lp(σ1, σ2; p=1) == 6.4
f_x(t::Real) = MarkovProcesses._f_step(x_obs, t_x, t)
......@@ -33,7 +33,7 @@ test_2 = isapprox(6.4, int; atol=err)
# Case 1 bis : inverse of case 1
test_1_bis = dist_lp(@view(y_obs[:]), @view(t_y[:]), @view(x_obs[:]), @view(t_x[:]); p=1) == 6.4
test_1_bis = dist_lp(@view(y_obs[:]), t_y, @view(x_obs[:]), t_x; p=1) == 6.4
test_1_bis = test_1_bis && dist_lp(σ2, σ1; p=1) == 6.4
return test_1 && test_1_bis && test_2
......
......@@ -14,8 +14,8 @@ for p = 1:2
d = dist_lp(σ1, σ2, "I"; p = p)
d2 = dist_lp(σ1, σ2; p = p)
f_x(t::Float64) = MarkovProcesses._f_step(σ1["I"], σ1["times"], t)
f_y(t::Float64) = MarkovProcesses._f_step(σ2["I"], σ2["times"], t)
f_x(t::Float64) = MarkovProcesses._f_step(σ1["I"], times(σ1), t)
f_y(t::Float64) = MarkovProcesses._f_step(σ2["I"], times(σ2), t)
diff_f(t) = abs(f_x(t) - f_y(t))^p
int_riemann = MarkovProcesses._riemann_sum(diff_f, SIR.t0, SIR.time_bound, 1E-3)
res_int_riemann = int_riemann^(1/p)
......
......@@ -6,7 +6,7 @@ load_model("ER")
σ = simulate(ER)
plt.figure()
plt.step(σ["times"], σ["P"], "ro--", marker="x", where="post", linewidth=1.0)
plt.step(times(σ), σ["P"], "ro--", marker="x", where="post", linewidth=1.0)
plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_er.png", dpi=480)
plt.close()
......
......@@ -8,7 +8,7 @@ ER_row_buffer.time_bound = 10.0
σ = _simulate_row_buffer(ER_row_buffer)
plt.figure()
plt.step(σ["times"], _get_values_row(σ,"P"), "ro--", marker="x", where="post", linewidth=1.0)
plt.step(times(σ), _get_values_row(σ,"P"), "ro--", marker="x", where="post", linewidth=1.0)
plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_er_row_buffer_bounded.png")
plt.close()
......
......@@ -6,7 +6,7 @@ load_model("SIR")
σ = simulate(SIR)
plt.figure()
plt.step(σ["times"], σ["I"], "ro--", marker="x", where="post", linewidth=1.0)
plt.step(times(σ), σ["I"], "ro--", marker="x", where="post", linewidth=1.0)
plt.savefig(get_module_path() * "/tests/simulation/res_pics/sim_sir.png")
plt.close()
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment