ParametricModel methods are implemented and tested.

Fix about methods of SynchronizedTrajectory + tests.

Simulation for ParametricModel works.

core/MarkovProcesses.jl

 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
 

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
 

core/lha.jl

 
-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]]
 

core/model.jl

 
 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.
 #

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
 

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(σ)

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")
 

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)
 

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
 

tests/simulation/sim_pm_er.jl

+
+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
+

tests/simulation/sim_pm_sync_er.jl

+
+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
+

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")

tests/unit/density_pm.jl

+
+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
+

tests/unit/draw_pm.jl

+
+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
+

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