From c0b8ba67d1be0ce504066bb6e902e236c48b28ec Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Fri, 27 Nov 2020 02:18:52 +0100
Subject: [PATCH] Reimplementation + methods of ParametricModel. It is more
clear and more efficient.
---
core/common.jl | 33 ++++++-----
core/model.jl | 117 +++++++++++++++++++++-----------------
tests/unit/model_prior.jl | 12 ++--
3 files changed, 91 insertions(+), 71 deletions(-)
diff --git a/core/common.jl b/core/common.jl
index a0026fb..a28d819 100644
--- a/core/common.jl
+++ b/core/common.jl
@@ -1,4 +1,7 @@
+import Distributions: Distribution, Univariate, Continuous, UnivariateDistribution,
+ MultivariateDistribution, product_distribution
+
abstract type Model end
abstract type AbstractTrajectory end
@@ -32,11 +35,6 @@ struct Trajectory <: AbstractTrajectory
transitions::Vector{Transition}
end
-struct ParametricModel
- m::Model
- map_l_param_dist::Dict{Vector{String},Distribution}
-end
-
struct Edge
transitions::Vector{Transition}
check_constraints::Function
@@ -76,6 +74,12 @@ struct SynchronizedTrajectory <: AbstractTrajectory
transitions::Vector{Transition}
end
+struct ParametricModel
+ m::Model
+ l_param::Vector{String}
+ dist::Distribution
+end
+
# Constructors
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,
@@ -105,16 +109,19 @@ 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 ParametricModel(m::Model, priors::Tuple{Vector{String},Distribution}...)
- map = Dict{Vector{String},Distribution}()
+function ParametricModel(am::Model, priors::Tuple{String,UnivariateDistribution}...)
+ m = get_proba_model(am)
+ l_param = String[]
+ l_dist = Distribution{Univariate,Continuous}[]
for prior in priors
check_vars = true
- for var in prior[1]
- check_vars = check_vars && var in keys(get_proba_model(m).map_param_idx)
- end
- @assert check_vars
- map[prior[1]] = prior[2]
+ str_p = prior[1]
+ dist = prior[2]
+ @assert str_p in keys(m.map_param_idx)
+ push!(l_param, str_p)
+ push!(l_dist, dist)
end
- return ParametricModel(m, map)
+ return ParametricModel(m, l_param, product_distribution(l_dist))
end
+
diff --git a/core/model.jl b/core/model.jl
index 359133d..7147553 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -1,4 +1,7 @@
+import Random: rand, rand!
+import Distributions: pdf
+
load_model(name_model::String) = include(get_module_path() * "/models/" * name_model * ".jl")
function _resize_trajectory!(values::Vector{Vector{Int}}, times::Vector{Float64},
@@ -29,7 +32,11 @@ end
Simulates a model. If `m::SynchronizedModel`, the simulation is synchronized with a
Linear Hybrid Automaton.
"""
-function simulate(m::ContinuousTimeModel)
+function simulate(m::ContinuousTimeModel; p::Union{Nothing,Vector{Float64}} = nothing)
+ p_sim = m.p
+ if p != nothing
+ p_sim = p
+ end
# First alloc
full_values = Vector{Vector{Int}}(undef, m.d)
for i = eachindex(full_values) full_values[i] = zeros(Int, m.estim_min_states) end
@@ -44,7 +51,7 @@ function simulate(m::ContinuousTimeModel)
xn = m.x0 # View for type stability
tn = m.t0
# at time n+1
- isabsorbing::Bool = m.isabsorbing(m.p,xn)
+ isabsorbing::Bool = m.isabsorbing(p_sim,xn)
# If x0 is absorbing
if isabsorbing
_resize_trajectory!(full_values, times, transitions, 1)
@@ -60,7 +67,7 @@ function simulate(m::ContinuousTimeModel)
l_tr = Transition[nothing]
# First we fill the allocated array
for i = 2:m.estim_min_states
- m.f!(vec_x, l_t, l_tr, xn, tn, m.p)
+ m.f!(vec_x, l_t, l_tr, xn, tn, p_sim)
tn = l_t[1]
if tn > m.time_bound
break
@@ -68,7 +75,7 @@ function simulate(m::ContinuousTimeModel)
n += 1
xn = vec_x
_update_values!(full_values, times, transitions, xn, tn, l_tr[1], i)
- isabsorbing = m.isabsorbing(m.p,xn)
+ isabsorbing = m.isabsorbing(p_sim,xn)
if isabsorbing
break
end
@@ -90,7 +97,7 @@ function simulate(m::ContinuousTimeModel)
i = 0
while i < m.buffer_size
i += 1
- m.f!(vec_x, l_t, l_tr, xn, tn, m.p)
+ m.f!(vec_x, l_t, l_tr, xn, tn, p_sim)
tn = l_t[1]
if tn > m.time_bound
i -= 1
@@ -99,7 +106,7 @@ function simulate(m::ContinuousTimeModel)
xn = vec_x
_update_values!(full_values, times, transitions,
xn, tn, l_tr[1], m.estim_min_states+size_tmp+i)
- isabsorbing = m.isabsorbing(m.p,xn)
+ isabsorbing = m.isabsorbing(p_sim,xn)
if isabsorbing
break
end
@@ -119,7 +126,11 @@ function simulate(m::ContinuousTimeModel)
end
return Trajectory(m, values, times, transitions)
end
-function simulate(product::SynchronizedModel)
+function simulate(product::SynchronizedModel; p::Union{Nothing,Vector{Float64}} = nothing)
+ p_sim = m.p
+ if p != nothing
+ p_sim = p
+ end
m = product.m
A = product.automaton
# First alloc
@@ -137,7 +148,7 @@ function simulate(product::SynchronizedModel)
xn = m.x0 # View for type stability
tn = m.t0
Sn = copy(S0)
- isabsorbing::Bool = m.isabsorbing(m.p,xn)
+ isabsorbing::Bool = m.isabsorbing(p_sim,xn)
isacceptedLHA::Bool = isaccepted(Sn)
# If x0 is absorbing
if isabsorbing || isacceptedLHA
@@ -155,7 +166,7 @@ function simulate(product::SynchronizedModel)
Snplus1 = copy(Sn)
# First we fill the allocated array
for i = 2:m.estim_min_states
- m.f!(vec_x, l_t, l_tr, xn, tn, m.p)
+ m.f!(vec_x, l_t, l_tr, xn, tn, p_sim)
tn = l_t[1]
if tn > m.time_bound
break
@@ -166,7 +177,7 @@ function simulate(product::SynchronizedModel)
next_state!(Snplus1, A, xn, tn, tr_n, Sn)
_update_values!(full_values, times, transitions, xn, tn, tr_n, i)
Sn = Snplus1
- isabsorbing = m.isabsorbing(m.p,xn)
+ isabsorbing = m.isabsorbing(p_sim,xn)
isacceptedLHA = isaccepted(Snplus1)
if isabsorbing || isacceptedLHA
break
@@ -189,7 +200,7 @@ function simulate(product::SynchronizedModel)
i = 0
while i < m.buffer_size
i += 1
- m.f!(vec_x, l_t, l_tr, xn, tn, m.p)
+ m.f!(vec_x, l_t, l_tr, xn, tn, p_sim)
tn = l_t[1]
if tn > m.time_bound
i -= 1
@@ -201,7 +212,7 @@ function simulate(product::SynchronizedModel)
_update_values!(full_values, times, transitions,
xn, tn, tr_n, m.estim_min_states+size_tmp+i)
Sn = Snplus1
- isabsorbing = m.isabsorbing(m.p,xn)
+ isabsorbing = m.isabsorbing(p_sim,xn)
isacceptedLHA = isaccepted(Snplus1)
if isabsorbing || isacceptedLHA
break
@@ -223,11 +234,24 @@ function simulate(product::SynchronizedModel)
return SynchronizedTrajectory(Sn, product, values, times, transitions)
end
-
function simulate(m::ContinuousTimeModel, n::Int)
end
-function set_observed_var!(m::Model, g::Vector{String})
+isbounded(m::ContinuousTimeModel) = m.time_bound < Inf
+function check_consistency(m::ContinuousTimeModel)
+ @assert m.d == length(m.map_var_idx)
+ @assert m.k == length(m.map_param_idx)
+ @assert m.k == length(m.p)
+ @assert length(m.g) <= m.d
+ @assert length(m._g_idx) == length(m.g)
+ @assert m.buffer_size >= 0
+ @assert typeof(m.isabsorbing(m.p, m.x0)) == Bool
+ return true
+end
+
+# Set and get Model fields
+function set_observed_var!(am::Model, g::Vector{String})
+ m = get_proba_model(am)
dobs = length(g)
_map_obs_var_idx = Dict{String}{Int}()
_g_idx = zeros(Int, dobs)
@@ -239,8 +263,8 @@ function set_observed_var!(m::Model, g::Vector{String})
m._g_idx = _g_idx
m._map_obs_var_idx = _map_obs_var_idx
end
-
-function observe_all!(m::Model)
+function observe_all!(am::Model)
+ m = get_proba_model(am)
g = Vector{String}(undef, m.d)
_g_idx = collect(1:m.d)
for var in keys(m.map_var_idx)
@@ -250,32 +274,17 @@ function observe_all!(m::Model)
m._g_idx = _g_idx
m._map_obs_var_idx = m.map_var_idx
end
-
-isbounded(m::ContinuousTimeModel) = m.time_bound < Inf
-function check_consistency(m::ContinuousTimeModel)
- @assert m.d == length(m.map_var_idx)
- @assert m.k == length(m.map_param_idx)
- @assert m.k == length(m.p)
- @assert length(m.g) <= m.d
- @assert length(m._g_idx) == length(m.g)
- @assert m.buffer_size >= 0
- @assert typeof(m.isabsorbing(m.p, m.x0)) == 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
+ for i = eachindex(l_name_p)
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)
+set_time_bound!(am::Model, b::Float64) = (get_proba_model(am).time_bound = b)
function getproperty(m::ContinuousTimeModel, sym::Symbol)
if sym == :dobs
@@ -284,33 +293,37 @@ function getproperty(m::ContinuousTimeModel, sym::Symbol)
return getfield(m, sym)
end
end
+
+function getproperty(pm::ParametricModel, sym::Symbol)
+ if sym == :df
+ return length(pm.l_param)
+ else
+ return getfield(pm, sym)
+ end
+end
get_proba_model(m::ContinuousTimeModel) = m
get_proba_model(sm::SynchronizedModel) = sm.m
-get_observed_var(m::ContinuousTimeModel) = m.g
-get_observed_var(sm::SynchronizedModel) = (sm.m).g
+get_proba_model(pm::ParametricModel) = get_proba_model(pm.m)
+get_observed_var(m::Model) = get_proba_model(am).g
# Prior methods
-function draw_model!(pm::ParametricModel)
- dict_dist = pm.map_l_param_dist
- for l_name in keys(dict_dist)
- if length(l_name) == 1
- set_param!(get_proba_model(pm.m), l_name[1], rand(dict_dist[l_name]))
- else
- set_param!(get_proba_model(pm.m), l_name, rand(dict_dist[l_name]))
- end
+draw_model!(pm::ParametricModel) = set_param!(get_proba_model(pm), pm.l_param, rand(pm.dist))
+draw!(vec_p::AbstractVector{Float64}, pm::ParametricModel) = rand!(pm.dist, vec_p)
+function draw!(mat_p::Matrix{Float64}, pm::ParametricModel)
+ for i = 1:size(mat_p)[2]
+ draw!(view(mat_p, :, i), pm)
end
end
-function dim_free_param(pm::ParametricModel)
- return 1
-end
-
-function draw!(vec_p::Vector{Float64}, pm::ParametricModel)
+prior_density(vec_p::AbstractVector{Float64}, pm::ParametricModel) = pdf(pm.dist, vec_p)
+function prior_density!(res_density::Vector{Float64}, mat_p::Matrix{Float64}, pm::ParametricModel)
+ for i = eachindex(res_density)
+ res_density[i] = prior_density(view(mat_p, :, i), pm)
+ end
end
-function draw!(mat_p::Matrix{Float64}, pm::ParametricModel)
-end
+function check_bounds() end
-function prior_density!(wl::Vector{Float64}, mat_p::Matrix{Float64}, pm::ParametricModel)
-end
+# to do: simulate(pm), create_res_dir, check_bounds, ajouter un champ _param_idx pour pm.
+#
diff --git a/tests/unit/model_prior.jl b/tests/unit/model_prior.jl
index b634aec..5d7721b 100644
--- a/tests/unit/model_prior.jl
+++ b/tests/unit/model_prior.jl
@@ -5,17 +5,17 @@ load_model("ER")
test_all = true
load_model("ER")
-prior1 = ParametricModel(ER, (["k2"], Uniform(2.0, 4.0)))
+prior1 = ParametricModel(ER, ("k2", Uniform(2.0, 4.0)))
draw_model!(prior1)
-test_all = test_all && 2.0 <= ER["k2"] <= 4.0
+test_all = test_all && 2.0 <= ER["k2"] <= 4.0 && prior1.df == 1
-prior2 = ParametricModel(ER, (["k3","k2"], Product(Uniform.([2.5,6.0], [3.5,7.0]))))
+prior2 = ParametricModel(ER, ["k3","k2"], Product(Uniform.([2.5,6.0], [3.5,7.0])))
draw_model!(prior2)
-test_all = test_all && 2.5 <= ER["k3"] <= 3.5 && 6.0 <= ER["k2"] <= 7.0
+test_all = test_all && 2.5 <= ER["k3"] <= 3.5 && 6.0 <= ER["k2"] <= 7.0 && prior2.df == 2
-prior3 = ParametricModel(ER, (["k3"], Uniform(10.0, 11.0)), (["k2"], Uniform(13.0, 14.0)))
+prior3 = ParametricModel(ER, ("k3", Uniform(10.0, 11.0)), ("k2", Uniform(13.0, 14.0)))
draw_model!(prior3)
-test_all = test_all && 10.0 <= ER["k3"] <= 11.0 && 13.0 <= ER["k2"] <= 14.0
+test_all = test_all && 10.0 <= ER["k3"] <= 11.0 && 13.0 <= ER["k2"] <= 14.0 && prior3.df == 2
return test_all
--
GitLab