From 938646900ddc77a3f4de323ff2321dc28b025cb6 Mon Sep 17 00:00:00 2001 From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr> Date: Fri, 27 Nov 2020 18:17:42 +0100 Subject: [PATCH] ParametricModel methods are implemented and tested. Fix about methods of SynchronizedTrajectory + tests. Simulation for ParametricModel works. --- core/MarkovProcesses.jl | 17 ++--- core/common.jl | 15 ++++- core/lha.jl | 2 - core/model.jl | 73 +++++++++++++++++----- core/plots.jl | 3 +- core/trajectory.jl | 12 ++-- core/utils.jl | 2 + tests/automaton_abc/R1.jl | 2 +- tests/run_simulation.jl | 2 + tests/simulation/sim_pm_er.jl | 15 +++++ tests/simulation/sim_pm_sync_er.jl | 17 +++++ tests/unit/check_trajectory_consistency.jl | 37 +++++++++-- tests/unit/density_pm.jl | 32 ++++++++++ tests/unit/draw_pm.jl | 24 +++++++ tests/unit/model_prior.jl | 18 +++--- 15 files changed, 220 insertions(+), 51 deletions(-) create mode 100644 tests/simulation/sim_pm_er.jl create mode 100644 tests/simulation/sim_pm_sync_er.jl create mode 100644 tests/unit/density_pm.jl create mode 100644 tests/unit/draw_pm.jl diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl index 914b006..0aa1edf 100644 --- a/core/MarkovProcesses.jl +++ b/core/MarkovProcesses.jl @@ -1,7 +1,8 @@ module MarkovProcesses import Base: +, -, * -import Base: copy, getfield, getindex, lastindex, setindex!, getproperty, setproperty! +import Base: copy, getfield, getindex, getproperty, lastindex +import Base: setindex!, setproperty!, fill! import StaticArrays: SVector import Distributions: Distribution, Product, Uniform, Normal @@ -21,27 +22,27 @@ export check_consistency, issteadystate, isaccepted # LHA related methods export init_state, next_state!, read_trajectory -export load_automaton, get_index, get_value, length_var, isaccepted +export get_index, get_value, length_var, isaccepted # Model related methods -export simulate, set_param!, get_param, set_observed_var!, observe_all! -export set_time_bound!, getproperty, draw!, draw_model! +export simulate, set_param!, set_time_bound!, set_observed_var!, observe_all! +export get_param, getproperty, get_proba_model, get_observed_var export isbounded, isaccepted, check_consistency -export load_model, get_module_path +export draw_model!, draw!, fill!, prior_pdf!, prior_pdf, insupport # Utils -export get_module_path, cosmos_get_values, load_plots +export get_module_path, cosmos_get_values +export load_model, load_automaton, load_plots # Algorithms export automaton_abc include("common.jl") - include("trajectory.jl") include("lha.jl") include("model.jl") include("utils.jl") -include(get_module_path() * "/algorithms/automaton_abc.jl") +include("../algorithms/automaton_abc.jl") end diff --git a/core/common.jl b/core/common.jl index a28d819..3c87419 100644 --- a/core/common.jl +++ b/core/common.jl @@ -78,6 +78,7 @@ struct ParametricModel m::Model l_param::Vector{String} dist::Distribution + _param_idx::Vector{Int} end # Constructors @@ -113,15 +114,25 @@ function ParametricModel(am::Model, priors::Tuple{String,UnivariateDistribution} m = get_proba_model(am) l_param = String[] l_dist = Distribution{Univariate,Continuous}[] + _param_idx = zeros(Int, 0) for prior in priors check_vars = true 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) + push!(l_dist, dist) + push!(_param_idx, m.map_param_idx[str_p]) end - return ParametricModel(m, l_param, product_distribution(l_dist)) + return ParametricModel(am, l_param, product_distribution(l_dist), _param_idx) end +function ParametricModel(am::Model, l_param::Vector{String}, dist::MultivariateDistribution) + m = get_proba_model(am) + _param_idx = zeros(Int, 0) + for str_p in l_param + push!(_param_idx, m.map_param_idx[str_p]) + end + return ParametricModel(am, l_param, dist, _param_idx) +end diff --git a/core/lha.jl b/core/lha.jl index 9e146da..e6f6680 100644 --- a/core/lha.jl +++ b/core/lha.jl @@ -1,6 +1,4 @@ -load_automaton(automaton::String) = include(get_module_path() * "/automata/$(automaton).jl") - length_var(A::LHA) = length(A.map_var_automaton_idx) get_value(A::LHA, x::Vector{Int}, var::String) = x[A.map_var_model_idx[var]] diff --git a/core/model.jl b/core/model.jl index 7147553..347716b 100644 --- a/core/model.jl +++ b/core/model.jl @@ -1,8 +1,6 @@ import Random: rand, rand! -import Distributions: pdf - -load_model(name_model::String) = include(get_module_path() * "/models/" * name_model * ".jl") +import Distributions: insupport, pdf function _resize_trajectory!(values::Vector{Vector{Int}}, times::Vector{Float64}, transitions::Vector{Transition}, size::Int) @@ -32,7 +30,7 @@ end Simulates a model. If `m::SynchronizedModel`, the simulation is synchronized with a Linear Hybrid Automaton. """ -function simulate(m::ContinuousTimeModel; p::Union{Nothing,Vector{Float64}} = nothing) +function simulate(m::ContinuousTimeModel; p::Union{Nothing,AbstractVector{Float64}} = nothing) p_sim = m.p if p != nothing p_sim = p @@ -126,13 +124,14 @@ function simulate(m::ContinuousTimeModel; p::Union{Nothing,Vector{Float64}} = no end return Trajectory(m, values, times, transitions) end -function simulate(product::SynchronizedModel; p::Union{Nothing,Vector{Float64}} = nothing) + +function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Float64}} = nothing) + m = product.m + A = product.automaton p_sim = m.p if p != nothing p_sim = p end - m = product.m - A = product.automaton # 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 @@ -234,10 +233,24 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,Vector{Float64}} return SynchronizedTrajectory(Sn, product, values, times, transitions) end +""" + `simulate(pm::ParametricModel, p_prior::AbstractVector{Float64}) + +Simulates the model contained in pm with p_prior values. +It simulates the model by taking the parameters contained in get_proba_model(pm).p and +replace the 1D parameters pm.l_param with p_prior. +""" +function simulate(pm::ParametricModel, p_prior::AbstractVector{Float64}) + full_p = copy(get_proba_model(pm).p) + full_p[pm._param_idx] = p_prior + + return simulate(pm.m; p = full_p) +end + function simulate(m::ContinuousTimeModel, n::Int) end -isbounded(m::ContinuousTimeModel) = m.time_bound < Inf +isbounded(m::Model) = get_proba_model(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) @@ -293,7 +306,6 @@ 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) @@ -301,28 +313,57 @@ function getproperty(pm::ParametricModel, sym::Symbol) return getfield(pm, sym) end end + get_proba_model(m::ContinuousTimeModel) = m get_proba_model(sm::SynchronizedModel) = sm.m get_proba_model(pm::ParametricModel) = get_proba_model(pm.m) + get_observed_var(m::Model) = get_proba_model(am).g # Prior methods +""" + `draw_model!(pm::ParametricModel)` + +Draw a parameter from the prior disitribution defined in `pm::ParametricModel` +and save it in the model contained in `pm`. +""" draw_model!(pm::ParametricModel) = set_param!(get_proba_model(pm), pm.l_param, rand(pm.dist)) +""" + `draw!(vec_p, pm)` + +Draw a parameter from the prior distribution defined in pm and stores it in vec_p. +""" draw!(vec_p::AbstractVector{Float64}, pm::ParametricModel) = rand!(pm.dist, vec_p) +""" + `draw!(mat_p, pm)` + +Draw `size(mat_p)[2]` (number of columns of mat_p) parameters from the prior distribution +defined in pm and stores it in mat_p. +""" function draw!(mat_p::Matrix{Float64}, pm::ParametricModel) for i = 1:size(mat_p)[2] draw!(view(mat_p, :, i), pm) end end - -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) +""" + `prior_pdf(p_prior, pm)` + +Computes the density of the prior distribution defined in pm on argument p_prior. +""" +prior_pdf(pm::ParametricModel, p_prior::AbstractVector{Float64}) = pdf(pm.dist, p_prior) +""" + `prior_pdf(vec_res, mat_p, pm)` + +Computes the density of the prior distribution defined in pm on each column +ov mat_p. Stores it in mat_p. (`length(vec_res) == size(mat_p)[2]`) +""" +function prior_pdf!(res_pdf::Vector{Float64}, pm::ParametricModel, mat_p::Matrix{Float64}) + for i = eachindex(res_pdf) + res_pdf[i] = prior_pdf(pm, view(mat_p, :, i)) end end - -function check_bounds() end +fill!(pm::ParametricModel, p_prior::Vector{Float64}) = get_proba_model(pm).p[pm._param_idx] = p_prior +insupport(pm::ParametricModel, p_prior::Vector{Float64}) = insupport(pm.dist, p_prior) # to do: simulate(pm), create_res_dir, check_bounds, ajouter un champ _param_idx pour pm. # diff --git a/core/plots.jl b/core/plots.jl index 4e50cfd..d174d15 100644 --- a/core/plots.jl +++ b/core/plots.jl @@ -22,7 +22,7 @@ function plot(σ::AbstractTrajectory, vars::String...; filename::String = "", pl end # Plots - p = plot(title = "Trajectory", palette = :lightrainbow) + p = plot(title = "Trajectory", palette = :lightrainbow, dpi = 480) for var in to_plot @assert var in get_obs_var(σ) plot!(p, times(σ), σ[var], @@ -47,7 +47,6 @@ function plot(σ::AbstractTrajectory, vars::String...; filename::String = "", pl display(p) else png(p, filename) - close(p) end end diff --git a/core/trajectory.jl b/core/trajectory.jl index 9ceb4b0..730340c 100644 --- a/core/trajectory.jl +++ b/core/trajectory.jl @@ -151,7 +151,7 @@ end function check_consistency(σ::AbstractTrajectory) test_length_var = true - for i = 1:(σ.m).dobs + for i = 1:get_proba_model(σ.m).dobs test_length_i = (length(σ.values[1]) == length(σ.values[i])) test_length_var = test_length_var && test_length_i end @@ -159,23 +159,23 @@ function check_consistency(σ::AbstractTrajectory) (length(σ.times) == length(σ.values[1])) && test_length_var end - @assert length_obs_var(σ) == σ.m.dobs + @assert length_obs_var(σ) == get_proba_model(σ.m).dobs return true end -issteadystate(σ::AbstractTrajectory) = (σ.m).isabsorbing((σ.m).p, view(reshape(σ[end], 1, m.d), 1, :)) +issteadystate(σ::AbstractTrajectory) = get_proba_model(σ.m).isabsorbing(get_proba_model(σ.m).p, view(reshape(σ[end], 1, m.d), 1, :)) # Properties of the trajectory length_states(σ::AbstractTrajectory) = length(σ.times) length_obs_var(σ::AbstractTrajectory) = length(σ.values) -get_obs_var(σ::AbstractTrajectory) = (σ.m).g +get_obs_var(σ::AbstractTrajectory) = get_proba_model(σ.m).g isbounded(σ::AbstractTrajectory) = σ.transitions[end] == nothing && length_states(σ) >= 2 isaccepted(σ::SynchronizedTrajectory) = isaccepted(σ.S) # Access to trajectory values -get_var_values(σ::AbstractTrajectory, var::String) = σ.values[(σ.m)._map_obs_var_idx[var]] +get_var_values(σ::AbstractTrajectory, var::String) = σ.values[get_proba_model(σ.m)._map_obs_var_idx[var]] get_state(σ::AbstractTrajectory, idx::Int) = [σ.values[i][idx] for i = 1:length(σ.values)] # /!\ Creates an array get_value(σ::AbstractTrajectory, var::String, idx::Int) = get_var_values(σ, var)[idx] -# Operation @ +# Operation σ@t function get_state_from_time(σ::AbstractTrajectory, t::Float64) @assert t >= 0.0 l_t = times(σ) diff --git a/core/utils.jl b/core/utils.jl index ff0d6be..648adee 100644 --- a/core/utils.jl +++ b/core/utils.jl @@ -19,5 +19,7 @@ function cosmos_get_values(name_file::String) return dict_values end +load_model(name_model::String) = include(get_module_path() * "/models/" * name_model * ".jl") +load_automaton(automaton::String) = include(get_module_path() * "/automata/$(automaton).jl") load_plots() = include(get_module_path() * "/core/plots.jl") diff --git a/tests/automaton_abc/R1.jl b/tests/automaton_abc/R1.jl index 1e3da5a..bd9536c 100644 --- a/tests/automaton_abc/R1.jl +++ b/tests/automaton_abc/R1.jl @@ -5,7 +5,7 @@ load_model("ER") load_automaton("automaton_F") A_F = create_automaton_F(ER, 50.0, 75.0, 0.025, 0.05, "P") sync_ER = A_F*ER -pm_sync_ER = ParametricModel(sync_ER, (["k3"], Uniform(0.0, 100.0))) +pm_sync_ER = ParametricModel(sync_ER, ("k3", Uniform(0.0, 100.0))) automaton_abc(pm_sync_ER) diff --git a/tests/run_simulation.jl b/tests/run_simulation.jl index 3bc43c9..5dee590 100644 --- a/tests/run_simulation.jl +++ b/tests/run_simulation.jl @@ -12,5 +12,7 @@ if !isdir(str_dir_pics) mkdir(str_dir_pics) end @test include("simulation/sim_sir_row_buffer_bounded.jl") @test include("simulation/sim_er.jl") @test include("simulation/sim_er_row_buffer_bounded.jl") + @test include("simulation/sim_pm_er.jl") + @test include("simulation/sim_pm_sync_er.jl") end diff --git a/tests/simulation/sim_pm_er.jl b/tests/simulation/sim_pm_er.jl new file mode 100644 index 0000000..c208806 --- /dev/null +++ b/tests/simulation/sim_pm_er.jl @@ -0,0 +1,15 @@ + +using MarkovProcesses + +load_plots() +load_model("ER") + +pm_ER = ParametricModel(ER, ("k1", Uniform(0.0,100.0)), ("k2", Uniform(0.0,100.0))) + +prior_p = [0.2, 40.0] + +σ = simulate(pm_ER, prior_p) +MarkovProcesses.plot(σ; filename = get_module_path() * "/tests/simulation/res_pics/sim_pm_er_long.png") + +return true + diff --git a/tests/simulation/sim_pm_sync_er.jl b/tests/simulation/sim_pm_sync_er.jl new file mode 100644 index 0000000..226dcb3 --- /dev/null +++ b/tests/simulation/sim_pm_sync_er.jl @@ -0,0 +1,17 @@ + +using MarkovProcesses + +load_plots() +load_model("ER") +load_automaton("automaton_F") + +A_F = create_automaton_F(ER, 0.0, 100.0, 7.0, 8.0, "P") +pm_sync_ER = ParametricModel(ER*A_F, ("k1", Uniform(0.0,100.0)), ("k2", Uniform(0.0,100.0))) + +prior_p = [0.2, 40.0] + +σ = simulate(pm_sync_ER, prior_p) +MarkovProcesses.plot(σ; filename = get_module_path() * "/tests/simulation/res_pics/sim_pm_sync_er_long.png") + +return true + diff --git a/tests/unit/check_trajectory_consistency.jl b/tests/unit/check_trajectory_consistency.jl index 1b3f39a..5ccecf5 100644 --- a/tests/unit/check_trajectory_consistency.jl +++ b/tests/unit/check_trajectory_consistency.jl @@ -3,13 +3,14 @@ using MarkovProcesses load_model("SIR") load_model("ER") +load_automaton("automaton_F") test_all = true nb_sim = 1000 function show_traj(io::IO, σ::AbstractTrajectory, m::Model) println(io, "length(σ.values[1]), length(σ.times), length(σ.transitions)") - println(io, "$(length(σ.values[1])), l$(length(σ.times)), $(length(σ.transitions))") + println(io, "$(length(σ.values[1])), $(length(σ.times)), $(length(σ.transitions))") println(io, "isbounded(m), isbounded(σ)") println(io, "$(isbounded(m)), $(isbounded(σ))") println(io, σ.values) @@ -30,15 +31,41 @@ for i = 1:nb_sim global test_all = test_all && check_consistency(σ) && check_consistency(σ2) && !isbounded(σ) && !isbounded(σ2) catch err - show_traj(stdout,σ, SIR) - show_traj(stdout,σ2, ER) + show_traj(stdout, σ, SIR) + show_traj(stdout, σ2, ER) + global σ_dump = σ + global σ2_dump = σ2 + throw(err) + end + if !test_all + show_traj(stdout, σ, SIR) + show_traj(stdout, σ2, ER) + global σ_dump = σ + global σ2_dump = σ2 + error("Ouch") + end +end + +new_SIR = deepcopy(SIR) +sync_SIR = new_SIR * create_automaton_F(new_SIR, 0.0, Inf, 100.0, 110.0, "I") +new_ER = deepcopy(ER) +sync_ER = new_ER * create_automaton_F(new_ER, 0.0, 100.0, 4.0, 5.0, "P") +for i = 1:nb_sim + local σ = simulate(sync_SIR) + local σ2 = simulate(sync_ER) + try + global test_all = test_all && check_consistency(σ) && check_consistency(σ2) && + !isbounded(σ) && !isbounded(σ2) + catch err + show_traj(stdout, σ, sync_SIR) + show_traj(stdout, σ2, sync_ER) global σ_dump = σ global σ2_dump = σ2 throw(err) end if !test_all - show_traj(stdout,σ, SIR) - show_traj(stdout,σ2, ER) + show_traj(stdout, σ, sync_SIR) + show_traj(stdout, σ2, sync_ER) global σ_dump = σ global σ2_dump = σ2 error("Ouch") diff --git a/tests/unit/density_pm.jl b/tests/unit/density_pm.jl new file mode 100644 index 0000000..bdabfae --- /dev/null +++ b/tests/unit/density_pm.jl @@ -0,0 +1,32 @@ + +using MarkovProcesses +import Distributions: rand, pdf, product_distribution +load_model("SIR") + +tol = 0.0 + +dist = Uniform(0.0, 1.0) +pm = ParametricModel(SIR, ("kr", dist)) +test1 = !insupport(pm, [2.0]) + +dist = Normal() +pm = ParametricModel(SIR, ("kr", dist)) +test2 = isapprox(prior_pdf(pm, [0.05]), pdf(dist, 0.05); atol = tol) + + +mat_u = [[rand()] for i = 1:10] +vec_u = [mat_u[i][1] for i = 1:10] +_prior_pdf(x::Vector{Float64}) = prior_pdf(pm, x) +test3 = isapprox(_prior_pdf.(mat_u), pdf.(dist, vec_u); atol = tol) + +dist, dist2 = Normal(), Normal(1.0, 2.0) +prod_dist = product_distribution([dist, dist2]) +pm = ParametricModel(SIR, ("ki", dist), ("kr", dist2)) +mat_u = rand(2,10) +vec_u = [mat_u[:,i] for i = 1:10] +vec_res = zeros(10) +prior_pdf!(vec_res, pm, mat_u) +test4 = isapprox(vec_res, [pdf(prod_dist, vec_u[i]) for i = 1:10]; atol = tol) + +return test1 && test2 && test3 && test4 + diff --git a/tests/unit/draw_pm.jl b/tests/unit/draw_pm.jl new file mode 100644 index 0000000..c4c8cbd --- /dev/null +++ b/tests/unit/draw_pm.jl @@ -0,0 +1,24 @@ + +using MarkovProcesses +load_model("ER") + +k1 = ER["k1"] +dist_mv_unif = Product(Uniform.([2.5,6.0], [3.5,7.0])) +pm = ParametricModel(ER, ["k3","k2"], dist_mv_unif) +draw_model!(pm) +test1 = 2.5 <= ER["k3"] <= 3.5 && 6.0 <= ER["k2"] <= 7.0 && pm.df == 2 + +p_drawn = copy(ER.p) +p_drawn_bis = copy((pm.m).p) +draw_model!(pm) +test2 = 2.5 <= ER["k3"] <= 3.5 && 6.0 <= ER["k2"] <= 7.0 && pm.df == 2 && p_drawn == p_drawn_bis && ER.p != p_drawn + +vec_p = zeros(2) +draw!(vec_p, pm) +test3 = 2.5 <= vec_p[1] <= 3.5 && 6.0 <= vec_p[2] <= 7.0 && pm.df == 2 + +fill!(pm, vec_p) +test4 = ER.p[pm._param_idx] == vec_p && ER.p[[3,2]] == vec_p && ER.p[1] == k1 + +return test1 && test2 && test3 && test4 + diff --git a/tests/unit/model_prior.jl b/tests/unit/model_prior.jl index 5d7721b..f5b4b9a 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))) -draw_model!(prior1) -test_all = test_all && 2.0 <= ER["k2"] <= 4.0 && prior1.df == 1 +pm1 = ParametricModel(ER, ("k2", Uniform(2.0, 4.0))) +draw_model!(pm1) +test_all = test_all && 2.0 <= ER["k2"] <= 4.0 && pm1.df == 1 -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 && prior2.df == 2 +pm2 = ParametricModel(ER, ["k3","k2"], Product(Uniform.([2.5,6.0], [3.5,7.0]))) +draw_model!(pm2) +test_all = test_all && 2.5 <= ER["k3"] <= 3.5 && 6.0 <= ER["k2"] <= 7.0 && pm2.df == 2 -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 && prior3.df == 2 +pm3 = ParametricModel(ER, ("k3", Uniform(10.0, 11.0)), ("k2", Uniform(13.0, 14.0))) +draw_model!(pm3) +test_all = test_all && 10.0 <= ER["k3"] <= 11.0 && 13.0 <= ER["k2"] <= 14.0 && pm3.df == 2 return test_all -- GitLab