ModelPrior => ParametricModel (more consistent)

First bricks of the automaton_abc algorithm.

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
 

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
 

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
+

core/plots.jl

 
 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)

core/utils.jl

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

tests/automaton_abc/R1.jl

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

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
 

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

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