From 61fa9d1c82f9894f691277959108e884cf9b5ba8 Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Thu, 3 Dec 2020 10:11:15 +0100
Subject: [PATCH] Reorganization of abc related algorithms. ABC SMC with
observations is implemented but not yet tested. automaton abc still works.
---
algorithms/{automaton_abc.jl => abc_smc.jl} | 66 +++++++++++++++++----
1 file changed, 54 insertions(+), 12 deletions(-)
rename algorithms/{automaton_abc.jl => abc_smc.jl} (71%)
diff --git a/algorithms/automaton_abc.jl b/algorithms/abc_smc.jl
similarity index 71%
rename from algorithms/automaton_abc.jl
rename to algorithms/abc_smc.jl
index c12cecb..adc2426 100644
--- a/algorithms/automaton_abc.jl
+++ b/algorithms/abc_smc.jl
@@ -41,17 +41,56 @@ function automaton_abc(pm::ParametricModel; nbr_particles::Int = 100, alpha::Flo
write(file_cfg, "kernel type : $(kernel_type) \n")
close(file_cfg)
if nprocs() == 1
- return _automaton_abc(pm, nbr_particles, alpha, kernel_type, NT, duration_time, bound_sim, dir_res, str_var_aut)
+ return _abc_smc(pm, nbr_particles, alpha, kernel_type, NT, duration_time, bound_sim, dir_res, str_var_aut)
end
- return _distributed_automaton_abc(pm, nbr_particles, alpha, kernel_type, NT, duration_time, bound_sim, dir_res, str_var_aut)
+ return _distributed_abc_smc(pm, nbr_particles, alpha, kernel_type, NT, duration_time, bound_sim, dir_res, str_var_aut)
end
+"""
+ `abc_smc(pm::ParametricModel, l_obs, func_dist; nbr_particles, alpha, kernel_type, NT
+ duration_tiùe, bound_sim, str_var_aut, verbose)`
+
+Run the ABC-SMC algorithm with the pm parametric model.
+`func_dist(l_sim, l_obs)` is the distance function between simulations and observation,
+it corresponds to \$\rho(\eta(y_sim), \eta(y_exp))\$.
+`l_obs::Vector{<:T2}` is a collection of observations.
+`dist` must have a signature of the form `func_dist(l_sim::Vector{T1}, l_obs::Vector{T2})`.
+If pm is defined on a ContinuousTimeModel, then `T1` should verify `T1 <: Trajectory`.
+
+!!! Distance function and distributed ABC
+
+ If you use `abc_smc` with multiple workers, `dist` has to be defined
+ on every workers by using @everywhere.
+
+"""
+function abc_smc(pm::ParametricModel, l_obs::AbstractVector, func_dist::Function;
+ nbr_particles::Int = 100, alpha::Float64 = 0.75, kernel_type = "mvnormal",
+ NT::Float64 = nbr_particles/2, duration_time::Float64 = Inf,
+ bound_sim::Int = typemax(Int), str_var_aut::String = "d", verbose::Int = 0)
+ @assert 0 < nbr_particles
+ @assert 0.0 < alpha < 1.0
+ @assert kernel_type in ["mvnormal", "knn_mvnormal"]
+ dir_res = create_results_dir()
+ file_cfg = open(dir_res * "cfg_abc_pmc.txt", "w")
+ write(file_cfg, "ParametricModel : $(pm) \n")
+ write(file_cfg, "Number of particles : $(nbr_particles) \n")
+ write(file_cfg, "alpha : $(alpha) \n")
+ write(file_cfg, "kernel type : $(kernel_type) \n")
+ close(file_cfg)
+ if nprocs() == 1
+ return _abc_smc(pm, nbr_particles, alpha, kernel_type, NT, duration_time, bound_sim, dir_res, str_var_aut; l_obs = l_obs, func_dist = func_dist)
+ end
+ return _distributed_abc_smc(pm, l_obs, dist, nbr_particles, alpha, kernel_type, NT, duration_time, bound_sim, dir_res, str_var_aut; l_obs = l_obs, func_dist = func_dist)
+end
+
+
# To code:
#Â Pkg related: draw!, prior_density!, create_results_dir
-function _automaton_abc(pm::ParametricModel, nbr_particles::Int, alpha::Float64,
- kernel_type::String, NT::Float64, duration_time::Float64, bound_sim::Int,
- dir_res::String, str_var_aut::String)
+function _abc_smc(pm::ParametricModel, nbr_particles::Int, alpha::Float64,
+ kernel_type::String, NT::Float64, duration_time::Float64,
+ bound_sim::Int, dir_res::String, str_var_aut::String;
+ l_obs::Union{Nothing,AbstractVector} = nothing, func_dist::Union{Nothing,Function} = nothing)
@info "ABC PMC with $(nworkers()) processus and $(Threads.nthreads()) threads"
begin_time = time()
nbr_p = pm.df
@@ -63,7 +102,7 @@ function _automaton_abc(pm::ParametricModel, nbr_particles::Int, alpha::Float64,
vec_dist = zeros(nbr_particles)
wl_old = zeros(nbr_particles)
@info "Step 1 : Init"
- _init_abc_automaton!(mat_p_old, vec_dist, pm, str_var_aut)
+ _init_abc_automaton!(mat_p_old, vec_dist, pm, str_var_aut; l_obs = l_obs, func_dist = func_dist)
prior_pdf!(wl_old, pm, mat_p_old)
normalize!(wl_old, 1)
ess = effective_sample_size(wl_old)
@@ -102,7 +141,8 @@ function _automaton_abc(pm::ParametricModel, nbr_particles::Int, alpha::Float64,
end
Threads.@threads for i = eachindex(vec_dist)
_update_param!(mat_p, vec_dist, l_nbr_sim, wl_current, i, pm, epsilon,
- wl_old, mat_p_old, mat_cov, tree_mat_p, kernel_type, str_var_aut)
+ wl_old, mat_p_old, mat_cov, tree_mat_p, kernel_type, str_var_aut;
+ l_obs = l_obs, func_dist = func_dist)
end
normalize!(wl_current, 1)
step_nbr_sim = sum(l_nbr_sim)
@@ -140,9 +180,10 @@ function _automaton_abc(pm::ParametricModel, nbr_particles::Int, alpha::Float64,
return r
end
-function _distributed_automaton_abc(pm::ParametricModel, nbr_particles::Int, alpha::Float64,
- kernel_type::String, NT::Float64, duration_time::Float64, bound_sim::Int,
- dir_res::String, str_var_aut::String)
+function _distributed_abc_smc(pm::ParametricModel, nbr_particles::Int, alpha::Float64,
+ kernel_type::String, NT::Float64, duration_time::Float64, bound_sim::Int,
+ dir_res::String, str_var_aut::String;
+ l_obs::Union{Nothing,AbstractVector} = nothing, func_dist::Union{Nothing,Function} = nothing)
@info "Distributed ABC PMC with $(nworkers()) processus and $(Threads.nthreads()) threads"
begin_time = time()
nbr_p = pm.df
@@ -154,7 +195,7 @@ function _distributed_automaton_abc(pm::ParametricModel, nbr_particles::Int, alp
vec_dist = zeros(nbr_particles)
wl_old = zeros(nbr_particles)
@info "Step 1 : Init"
- _init_abc_automaton!(mat_p_old, vec_dist, pm, str_var_aut)
+ _init_abc_automaton!(mat_p_old, vec_dist, pm, str_var_aut; l_obs = l_obs, func_dist = func_dist)
prior_pdf!(wl_old, pm, mat_p_old)
normalize!(wl_old, 1)
ess = effective_sample_size(wl_old)
@@ -199,7 +240,8 @@ function _distributed_automaton_abc(pm::ParametricModel, nbr_particles::Int, alp
@async l_nbr_sim[t_id_w] =
remotecall_fetch(() -> _task_worker!(d_mat_p, d_vec_dist, d_wl_current,
pm, epsilon, wl_old, mat_p_old, mat_cov, tree_mat_p,
- kernel_type, str_var_aut), id_w)
+ kernel_type, str_var_aut;
+ l_obs = l_obs, func_dist = func_dist), id_w)
end
wl_current = convert(Array, d_wl_current)
normalize!(wl_current, 1)
--
GitLab