diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl index e547640d4d41a7a0b2c7645e08e8786cf11d9074..914b006e458f6b83ce5bfb87797c377a8ad5dd45 100644 --- a/core/MarkovProcesses.jl +++ b/core/MarkovProcesses.jl @@ -10,7 +10,7 @@ export Distribution, Product, Uniform, Normal # Common types and constructors export Observations, AbstractTrajectory, Trajectory, SynchronizedTrajectory -export Model, ContinuousTimeModel, SynchronizedModel, ModelPrior +export Model, ContinuousTimeModel, SynchronizedModel, ParametricModel export LHA, StateLHA, Edge # Trajectory related methods @@ -25,19 +25,23 @@ export load_automaton, 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! +export set_time_bound!, getproperty, draw!, draw_model! export isbounded, isaccepted, check_consistency export load_model, get_module_path # Utils export get_module_path, cosmos_get_values, 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") end diff --git a/core/common.jl b/core/common.jl index be1a32f8eadffdc9a01fa99a50277106f881e381..a0026fb5fbc84b183e90fbb77fd8b1a630e55c67 100644 --- a/core/common.jl +++ b/core/common.jl @@ -32,7 +32,7 @@ struct Trajectory <: AbstractTrajectory transitions::Vector{Transition} end -struct ModelPrior +struct ParametricModel m::Model map_l_param_dist::Dict{Vector{String},Distribution} end @@ -105,11 +105,16 @@ 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}...) +function ParametricModel(m::Model, priors::Tuple{Vector{String},Distribution}...) map = Dict{Vector{String},Distribution}() 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] end - return ModelPrior(m, map) + return ParametricModel(m, map) end diff --git a/core/model.jl b/core/model.jl index 0a305c3fcb878f4fb2ddae6add7d099190ed1e8a..359133dba01e47cf777da21dcc55888716da7efc 100644 --- a/core/model.jl +++ b/core/model.jl @@ -286,16 +286,31 @@ function getproperty(m::ContinuousTimeModel, sym::Symbol) 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 # Prior methods -function draw!(Π::ModelPrior) - dict_dist = Π.map_l_param_dist +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(Π.m), l_name[1], rand(dict_dist[l_name])) + set_param!(get_proba_model(pm.m), l_name[1], rand(dict_dist[l_name])) else - set_param!(get_proba_model(Π.m), l_name, rand(dict_dist[l_name])) + set_param!(get_proba_model(pm.m), l_name, rand(dict_dist[l_name])) end end end +function dim_free_param(pm::ParametricModel) + return 1 +end + +function draw!(vec_p::Vector{Float64}, pm::ParametricModel) +end + +function draw!(mat_p::Matrix{Float64}, pm::ParametricModel) +end + +function prior_density!(wl::Vector{Float64}, mat_p::Matrix{Float64}, pm::ParametricModel) +end + diff --git a/core/plots.jl b/core/plots.jl index b2092ef7f4468ef042a475dbc13e8369ed2d3341..4e50cfd2444d9d96dd7999f0dcf9056a02e0b6e2 100644 --- a/core/plots.jl +++ b/core/plots.jl @@ -1,7 +1,16 @@ import Plots: plot, plot!, scatter! -import Plots: palette, display, png, close +import Plots: current, palette, display, png, close +""" + `plot(σ, var...; plot_transitions=false)` + +Plot a simulated trajectory σ. var... is a tuple of stirng variables. +`plot(σ)` will plot all the variables simulated in σ +whereas `plot(σ, "I", "R")` only plots the variables I and R of the trajectory (if it exists). +If `plot_transitions=true`, a marker that corresponds to a transition of the model will be plotted +at each break of the trajectory. +""" function plot(σ::AbstractTrajectory, vars::String...; filename::String = "", plot_transitions = false) # Setup palette_tr = palette(:default) diff --git a/core/utils.jl b/core/utils.jl index c403635c85567a464600399d7d9eba0875c4aa6a..ff0d6bea939575a8db5084a9ace3383164e1f307 100644 --- a/core/utils.jl +++ b/core/utils.jl @@ -1,6 +1,10 @@ get_module_path() = dirname(dirname(pathof(@__MODULE__))) +function create_results_dir() + return "./" +end + function cosmos_get_values(name_file::String) output_file = open(name_file) dict_values = Dict{String}{Float64}() diff --git a/tests/automaton_abc/R1.jl b/tests/automaton_abc/R1.jl new file mode 100644 index 0000000000000000000000000000000000000000..1e3da5aa54d89ffe11f43555903b4dc802a4cfa9 --- /dev/null +++ b/tests/automaton_abc/R1.jl @@ -0,0 +1,11 @@ + +using MarkovProcesses + +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))) + +automaton_abc(pm_sync_ER) + diff --git a/tests/run_unit.jl b/tests/run_unit.jl index a384dcb288e0879e736960ecead45c40be5dee84..e67bfff5bc391c1ef70d380a5dc7f51832b693ab 100644 --- a/tests/run_unit.jl +++ b/tests/run_unit.jl @@ -2,27 +2,36 @@ using Test @testset "Unit tests" begin - @test include("unit/load_model.jl") - @test include("unit/load_model_bench.jl") - @test include("unit/load_module.jl") - @test include("unit/simulate_sir.jl") - @test include("unit/simulate_sir_bounded.jl") - @test include("unit/simulate_er.jl") + @test include("unit/absorbing_x0.jl") + @test include("unit/change_obs_var_sir.jl") @test include("unit/change_obs_var_sir_2.jl") + @test include("unit/check_model_consistency.jl") + @test include("unit/check_trajectory_consistency.jl") + @test include("unit/create_automata.jl") + @test include("unit/create_models.jl") + + @test include("unit/dist_lp.jl") + @test include("unit/dist_lp_var.jl") + @test include("unit/getindex_access_trajectory.jl") @test include("unit/is_always_bounded_sir.jl") - @test include("unit/length_obs_var.jl") - @test include("unit/dist_lp_var.jl") - @test include("unit/dist_lp.jl") + @test include("unit/l_dist_lp.jl") - @test include("unit/check_trajectory_consistency.jl") - @test include("unit/check_model_consistency.jl") + @test include("unit/length_obs_var.jl") + @test include("unit/load_model.jl") + @test include("unit/load_model_bench.jl") + @test include("unit/load_module.jl") + @test include("unit/long_sim_er.jl") + + @test include("unit/model_prior.jl") + @test include("unit/observe_all.jl") + + @test include("unit/set_param.jl") @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") - @test include("unit/absorbing_x0.jl") + @test include("unit/simulate_sir.jl") + @test include("unit/simulate_sir_bounded.jl") + @test include("unit/simulate_er.jl") end diff --git a/tests/unit/model_prior.jl b/tests/unit/model_prior.jl index 0b04b511c521794b41397c00f92c246421931e1b..b634aec8df83e6c377b0ac546254fb91c8b5e8ee 100644 --- a/tests/unit/model_prior.jl +++ b/tests/unit/model_prior.jl @@ -5,16 +5,16 @@ load_model("ER") test_all = true load_model("ER") -prior1 = ModelPrior(ER, (["k2"], Uniform(2.0, 4.0))) -draw!(prior1) +prior1 = ParametricModel(ER, (["k2"], Uniform(2.0, 4.0))) +draw_model!(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) +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 -prior3 = ModelPrior(ER, (["k3"], Uniform(10.0, 11.0)), (["k2"], Uniform(13.0, 14.0))) -draw!(prior3) +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 return test_all diff --git a/tests/unit/observe_all.jl b/tests/unit/observe_all.jl index 3608bfc4ca812d598b81c92992b5ae15422afac4..ff3b5f3ca4fa37b40ff20b9f7ad9188108b2a045 100644 --- a/tests/unit/observe_all.jl +++ b/tests/unit/observe_all.jl @@ -7,5 +7,5 @@ load_model("SIR") observe_all!(ER) observe_all!(SIR) -return ER.g == ["E", "S", "ER", "P"] && SIR.g == ["S", "I", "R"] +return (ER.g == ["E", "S", "ES", "P"] && SIR.g == ["S", "I", "R"])