From d745815538d8d60a822789989ee9d24a91b4bc3b Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Tue, 24 Nov 2020 10:32:06 +0100
Subject: [PATCH] Implement of ModelPrior + test. Add several convenient
methods for parameters of model.
Next step: implementation of simulate with resize!.
---
bench/pygmalion/multiple_long_sim_er.jl | 9 ++-
core/MarkovProcesses.jl | 7 +-
core/common.jl | 18 ++++-
core/model.jl | 103 +++++++++++++++++-------
tests/run_unit.jl | 1 +
tests/unit/model_prior.jl | 21 +++++
tests/unit/set_param.jl | 20 +++++
7 files changed, 143 insertions(+), 36 deletions(-)
create mode 100644 tests/unit/model_prior.jl
diff --git a/bench/pygmalion/multiple_long_sim_er.jl b/bench/pygmalion/multiple_long_sim_er.jl
index 23244dc..18552e0 100644
--- a/bench/pygmalion/multiple_long_sim_er.jl
+++ b/bench/pygmalion/multiple_long_sim_er.jl
@@ -13,7 +13,7 @@ pygmalion.load_model(str_m)
str_oml = "P,R,time"
ll_om = split(str_oml, ",")
x0 = State(95.0, 95.0, 0.0, 0.0, 0.0, 0.0)
-p_true = Parameters(1.0, 1.0, 1.0)
+p_true = Parameters(0.2, 40.0, 1.0)
u = Control(20.0)
tml = 1:400
g_all = create_observation_function([ObserverModel(str_oml, tml)])
@@ -29,6 +29,9 @@ println("MarkovProcesses:")
using MarkovProcesses
MarkovProcesses.load_model("ER")
ER.time_bound = 20.0
+ER.estim_min_states = 7000
+set_param!(ER, "k1", 0.2)
+set_param!(ER, "k2", 40.0)
@timev MarkovProcesses.simulate(ER)
println("Default buffer size=10")
@@ -40,6 +43,6 @@ println("Buffer size 100")
ER.buffer_size = 100
b1_100 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate($ER) end
b2_100 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end
-@show minimum(b1), mean(b1), maximum(b1)
-@show minimum(b2), mean(b2), maximum(b2)
+@show minimum(b1_100), mean(b1_100), maximum(b1_100)
+@show minimum(b2_100), mean(b2_100), maximum(b2_100)
diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index 85eb0f6..11c2b2d 100644
--- a/core/MarkovProcesses.jl
+++ b/core/MarkovProcesses.jl
@@ -4,10 +4,13 @@ import Base: +, -, *
import Base: copy, getfield, getindex, lastindex, setindex!, getproperty, setproperty!
import StaticArrays: SVector
+import Distributions: Distribution, Product, Uniform, Normal
+
+export Distribution, Product, Uniform, Normal
# Common types and constructors
export Observations, AbstractTrajectory, Trajectory, SynchronizedTrajectory
-export Model, ContinuousTimeModel, SynchronizedModel
+export Model, ContinuousTimeModel, SynchronizedModel, ModelPrior
export LHA, StateLHA, Edge
# Trajectory related methods
@@ -22,7 +25,7 @@ export load_automaton, get_index, get_value, length_var, isaccepted
# Model related methods
export simulate, set_param!, get_param, set_observed_var!
-export set_time_bound!, getproperty
+export set_time_bound!, getproperty, draw!
export isbounded, isaccepted, check_consistency
export load_model, get_module_path
diff --git a/core/common.jl b/core/common.jl
index 8f7b015..be1a32f 100644
--- a/core/common.jl
+++ b/core/common.jl
@@ -22,6 +22,7 @@ mutable struct ContinuousTimeModel <: Model
isabsorbing::Function
time_bound::Float64
buffer_size::Int
+ estim_min_states::Int
end
struct Trajectory <: AbstractTrajectory
@@ -31,6 +32,11 @@ struct Trajectory <: AbstractTrajectory
transitions::Vector{Transition}
end
+struct ModelPrior
+ m::Model
+ map_l_param_dist::Dict{Vector{String},Distribution}
+end
+
struct Edge
transitions::Vector{Transition}
check_constraints::Function
@@ -74,7 +80,7 @@ end
function ContinuousTimeModel(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_transitions::Vector{String},
p::Vector{Float64}, x0::Vector{Int}, t0::Float64,
f!::Function, isabsorbing::Function;
- g::Vector{String} = keys(map_var_idx), time_bound::Float64 = Inf, buffer_size::Int = 10)
+ g::Vector{String} = keys(map_var_idx), time_bound::Float64 = Inf, buffer_size::Int = 10, estim_min_states::Int = 50)
dobs = length(g)
_map_obs_var_idx = Dict()
_g_idx = Vector{Int}(undef, dobs)
@@ -90,7 +96,7 @@ function ContinuousTimeModel(d::Int, k::Int, map_var_idx::Dict, map_param_idx::D
@warn "You have possibly redefined a function Model.isabsorbing used in a previously instantiated model."
end
- return ContinuousTimeModel(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_transitions, p, x0, t0, f!, g, _g_idx, isabsorbing, time_bound, buffer_size)
+ return ContinuousTimeModel(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_transitions, p, x0, t0, f!, g, _g_idx, isabsorbing, time_bound, buffer_size, estim_min_states)
end
LHA(A::LHA, map_var::Dict{String,Int}) = LHA(A.l_transitions, A.l_loc, A.Λ,
@@ -99,3 +105,11 @@ LHA(A::LHA, map_var::Dict{String,Int}) = LHA(A.l_transitions, A.l_loc, A.Λ,
Base.:*(m::ContinuousTimeModel, A::LHA) = SynchronizedModel(m, A)
Base.:*(A::LHA, m::ContinuousTimeModel) = SynchronizedModel(m, A)
+function ModelPrior(m::Model, priors::Tuple{Vector{String},Distribution}...)
+ map = Dict{Vector{String},Distribution}()
+ for prior in priors
+ map[prior[1]] = prior[2]
+ end
+ return ModelPrior(m, map)
+end
+
diff --git a/core/model.jl b/core/model.jl
index bbbf51d..e449830 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -4,10 +4,11 @@ load_model(name_model::String) = include(get_module_path() * "/models/" * name_m
# Simulation
function simulate(m::ContinuousTimeModel)
# trajectory fields
- full_values = Matrix{Int}(undef, 1, m.d)
- full_values[1,:] = m.x0
+ full_values = Vector{Vector{Int}}(undef, m.d)
+ for i = 1:m.d full_values[i] = Int[m.x0[i]] end
+ for i = 1:m.d sizehint!(full_values[i], m.estim_min_states) end
times = Float64[m.t0]
- transitions = Union{String,Nothing}[nothing]
+ transitions = Transition[nothing]
# values at time n
n = 0
xn = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
@@ -17,39 +18,47 @@ function simulate(m::ContinuousTimeModel)
l_t = zeros(Float64, m.buffer_size)
l_tr = Vector{String}(undef, m.buffer_size)
isabsorbing::Bool = m.isabsorbing(m.p,xn)
+ end_idx = -1
+ # use sizehint! ?
while !isabsorbing && tn <= m.time_bound
- i = 0
- while i < m.buffer_size && !isabsorbing && tn <= m.time_bound
- i += 1
+ for i = 1:m.buffer_size
m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
tn = l_t[i]
if tn > m.time_bound
i -= 1 # 0 is an ok value, 1:0 is allowed
+ end_idx = i
break
end
xn = view(mat_x, i, :)
isabsorbing = m.isabsorbing(m.p,xn)
+ if isabsorbing
+ end_idx = i
+ break
+ end
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))
- n += i
- #isabsorbing = m.isabsorbing(m.p,xn)
+ if end_idx != -1
+ break
+ end
+ for k = 1:m.d
+ append!(full_values[k], view(mat_x, :, k))
+ end
+ append!(times, l_t)
+ append!(transitions, l_tr)
+ n += m.buffer_size
+ end
+ for k = 1:m.d
+ append!(full_values[k], view(mat_x, 1:end_idx, k))
end
+ append!(times, view(l_t, 1:end_idx))
+ append!(transitions, view(l_tr, 1:end_idx))
if isbounded(m)
- #=
- if times[end] > m.time_bound
- full_values[end,:] = full_values[end-1,:]
- times[end] = m.time_bound
- transitions[end] = nothing
- else
+ for k = 1:m.d
+ push!(full_values[k], full_values[k][end])
end
- =#
- full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d))
push!(times, m.time_bound)
push!(transitions, nothing)
end
- values = view(full_values, :, m._g_idx)
+ values = full_values[m._g_idx]
return Trajectory(m, values, times, transitions)
end
@@ -57,12 +66,13 @@ function simulate(product::SynchronizedModel)
# trajectory fields
m = product.m
A = product.automaton
- full_values = Matrix{Int}(undef, 1, m.d)
- full_values[1,:] = m.x0
+ full_values = Vector{Vector{Int}}(undef, m.d)
+ for i = 1:m.d full_values[i] = Int[m.x0[i]] end
+ for i = 1:m.d sizehint!(full_values[i], m.estim_min_states) end
times = Float64[m.t0]
transitions = Union{String,Nothing}[nothing]
reshaped_x0 = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
- S0 = init_state(A, reshaped_x0, m.t0)
+ S0 = init_state(A, m.x0, m.t0)
# values at time n
n = 0
xn = reshaped_x0
@@ -75,7 +85,6 @@ function simulate(product::SynchronizedModel)
isabsorbing::Bool = m.isabsorbing(m.p,xn)
isacceptedLHA::Bool = isaccepted(Sn)
Snplus1 = copy(Sn)
- Sn_dump = Sn
while !isabsorbing && tn <= m.time_bound && !isacceptedLHA
i = 0
while i < m.buffer_size && !isabsorbing && tn <= m.time_bound && !isacceptedLHA
@@ -93,7 +102,9 @@ function simulate(product::SynchronizedModel)
isabsorbing = m.isabsorbing(m.p,xn)
isacceptedLHA = isaccepted(Snplus1)
end
- full_values = vcat(full_values, view(mat_x, 1:i, :))
+ for k = 1:m.d
+ append!(full_values[k], view(mat_x, 1:i, k))
+ end
append!(times, view(l_t, 1:i))
append!(transitions, view(l_tr, 1:i))
n += i
@@ -101,11 +112,13 @@ function simulate(product::SynchronizedModel)
# When the trajectory is accepted, we should not add an end value
if isbounded(m) && !isaccepted(Sn)
@assert times[end] < m.time_bound
- #full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d))
- #push!(times, m.time_bound)
- #push!(transitions, nothing)
+ for k = 1:m.d
+ push!(full_values[k], full_values[k][end])
+ end
+ push!(times, m.time_bound)
+ push!(transitions, nothing)
end
- values = view(full_values, :, m._g_idx)
+ values = full_values[m._g_idx]
return SynchronizedTrajectory(Snplus1, product, values, times, transitions)
end
@@ -136,8 +149,40 @@ function check_consistency(m::ContinuousTimeModel)
@assert typeof(m.isabsorbing(m.p, view(reshape(m.x0, 1, m.d), 1, :))) == Bool
return true
end
+
set_param!(m::ContinuousTimeModel, p::Vector{Float64}) = (m.p = p)
set_param!(m::ContinuousTimeModel, name_p::String, p_i::Float64) = (m.p[m.map_param_idx[name_p]] = p_i)
+function set_param!(m::ContinuousTimeModel, l_name_p::Vector{String}, p::Vector{Float64})
+ nb_param = length(l_name_p)
+ for i = 1:nb_param
+ set_param!(m, l_name_p[i], p[i])
+ end
+end
+
get_param(m::ContinuousTimeModel) = m.p
+getindex(m::ContinuousTimeModel, name_p::String) = m.p[m.map_param_idx[name_p]]
set_time_bound!(m::ContinuousTimeModel, b::Float64) = (m.time_bound = b)
+set_time_bound!(sm::SynchronizedModel, b::Float64) = set_time_bound!(sm.m, b)
+
+function getproperty(m::ContinuousTimeModel, sym::Symbol)
+ if sym == :dobs
+ return length(m.g)
+ else
+ return getfield(m, sym)
+ end
+end
+get_proba_model(m::ContinuousTimeModel) = m
+get_proba_model(sm::SynchronizedModel) = sm.m
+
+# Prior methods
+function draw!(Π::ModelPrior)
+ dict_dist = Π.map_l_param_dist
+ for l_name in keys(dict_dist)
+ if length(l_name) == 1
+ set_param!(get_proba_model(Π.m), l_name[1], rand(dict_dist[l_name]))
+ else
+ set_param!(get_proba_model(Π.m), l_name, rand(dict_dist[l_name]))
+ end
+ end
+end
diff --git a/tests/run_unit.jl b/tests/run_unit.jl
index d81c2e7..3202f53 100644
--- a/tests/run_unit.jl
+++ b/tests/run_unit.jl
@@ -22,5 +22,6 @@ using Test
@test include("unit/side_effects_1.jl")
@test include("unit/create_models.jl")
@test include("unit/create_automata.jl")
+ @test include("unit/model_prior.jl")
end
diff --git a/tests/unit/model_prior.jl b/tests/unit/model_prior.jl
new file mode 100644
index 0000000..0b04b51
--- /dev/null
+++ b/tests/unit/model_prior.jl
@@ -0,0 +1,21 @@
+
+using MarkovProcesses
+load_model("ER")
+
+test_all = true
+
+load_model("ER")
+prior1 = ModelPrior(ER, (["k2"], Uniform(2.0, 4.0)))
+draw!(prior1)
+test_all = test_all && 2.0 <= ER["k2"] <= 4.0
+
+prior2 = ModelPrior(ER, (["k3","k2"], Product(Uniform.([2.5,6.0], [3.5,7.0]))))
+draw!(prior2)
+test_all = test_all && 2.5 <= ER["k3"] <= 3.5 && 6.0 <= ER["k2"] <= 7.0
+
+prior3 = ModelPrior(ER, (["k3"], Uniform(10.0, 11.0)), (["k2"], Uniform(13.0, 14.0)))
+draw!(prior3)
+test_all = test_all && 10.0 <= ER["k3"] <= 11.0 && 13.0 <= ER["k2"] <= 14.0
+
+return test_all
+
diff --git a/tests/unit/set_param.jl b/tests/unit/set_param.jl
index 222c845..f99c2f9 100644
--- a/tests/unit/set_param.jl
+++ b/tests/unit/set_param.jl
@@ -22,6 +22,11 @@ k3 = 10.0
set_param!(ER, "k3", 10.0)
test_all = test_all && (ER.p == [0.5, 4.0, 10.0])
+l_k = [20.0, 10.0]
+l_name = ["k3", "k2"]
+set_param!(ER, l_name, l_k)
+test_all = test_all && (ER.p == [0.5, 10.0, 20.0])
+
new_param_SIR = [0.0013, 0.08]
set_param!(SIR, new_param_SIR)
test_all = test_all && (SIR.p == new_param_SIR)
@@ -34,5 +39,20 @@ ki = 0.011
set_param!(SIR, "ki", ki)
test_all = test_all && (SIR.p == [0.011, 0.06])
+l_name = ["kr", "ki"]
+l_p = [0.03, 0.015]
+set_param!(SIR, l_name, l_p)
+test_all = test_all && (SIR.p == [0.015, 0.03])
+
+l_name = ["ki", "kr"]
+l_p = [0.03, 0.015]
+set_param!(SIR, l_name, l_p)
+test_all = test_all && (SIR.p == [0.03, 0.015])
+
+l_name = ["kr"]
+l_p = [0.02]
+set_param!(SIR, l_name, l_p)
+test_all = test_all && (SIR.p == [0.03, 0.02])
+
return test_all
--
GitLab