Reorganization of abc related algorithms. ABC SMC with observations is...

Reorganization of abc related algorithms. ABC SMC with observations is implemented but not yet tested. automaton abc still works.

algorithms/automaton_abc.jl → 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)