Commit 4c4abbfa authored by Bentriou Mahmoud's avatar Bentriou Mahmoud
Browse files

Small fix in ABC euclideau distance automaton for distributed abc.

Fix of distributed abc smc tests + add of a test.
Preparation for the new feature: add a stop criteria in the simulation.
All tests passed.
parent de9fbd7e
...@@ -125,7 +125,6 @@ function _abc_smc(pm::ParametricModel, nbr_particles::Int, tolerance::Float64, a ...@@ -125,7 +125,6 @@ function _abc_smc(pm::ParametricModel, nbr_particles::Int, tolerance::Float64, a
l_nbr_sim = zeros(Int, nbr_particles) l_nbr_sim = zeros(Int, nbr_particles)
while (epsilon > last_epsilon) && (current_time - begin_time <= duration_time) && (nbr_tot_sim <= bound_sim) while (epsilon > last_epsilon) && (current_time - begin_time <= duration_time) && (nbr_tot_sim <= bound_sim)
t += 1 t += 1
@show pm.m.automaton.ϵ, epsilon
begin_time_ite = time() begin_time_ite = time()
@info "Step $t" @info "Step $t"
# Set new epsilon # Set new epsilon
......
...@@ -25,8 +25,8 @@ end ...@@ -25,8 +25,8 @@ end
end end
# A push! method implementend by myself because of preallocation of edge_candidates # A push! method implementend by myself because of preallocation of edge_candidates
function _push_edge!(edge_candidates::Vector{<:EdgeABCEuclideanDistanceAutomaton}, @everywhere function _push_edge!(edge_candidates::Vector{<:EdgeABCEuclideanDistanceAutomaton},
edge::EdgeABCEuclideanDistanceAutomaton, nbr_candidates::Int) edge::EdgeABCEuclideanDistanceAutomaton, nbr_candidates::Int)
if nbr_candidates < length(edge_candidates) if nbr_candidates < length(edge_candidates)
edge_candidates[nbr_candidates+1] = edge edge_candidates[nbr_candidates+1] = edge
else else
...@@ -34,12 +34,12 @@ function _push_edge!(edge_candidates::Vector{<:EdgeABCEuclideanDistanceAutomaton ...@@ -34,12 +34,12 @@ function _push_edge!(edge_candidates::Vector{<:EdgeABCEuclideanDistanceAutomaton
end end
end end
function _find_edge_candidates!(edge_candidates::Vector{EdgeABCEuclideanDistanceAutomaton}, @everywhere function _find_edge_candidates!(edge_candidates::Vector{EdgeABCEuclideanDistanceAutomaton},
edges_from_current_loc::Dict{Location,Vector{EdgeABCEuclideanDistanceAutomaton}}, edges_from_current_loc::Dict{Location,Vector{EdgeABCEuclideanDistanceAutomaton}},
Λ::Dict{Location,InvariantPredicateFunction}, Λ::Dict{Location,InvariantPredicateFunction},
S_time::Float64, S_values::Vector{Float64}, S_time::Float64, S_values::Vector{Float64},
x::Vector{Int}, p::Vector{Float64}, ϵ::Float64, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64,
only_asynchronous::Bool) only_asynchronous::Bool)
nbr_candidates = 0 nbr_candidates = 0
for target_loc in keys(edges_from_current_loc) for target_loc in keys(edges_from_current_loc)
if !Λ[target_loc](x) continue end if !Λ[target_loc](x) continue end
...@@ -61,8 +61,8 @@ function _find_edge_candidates!(edge_candidates::Vector{EdgeABCEuclideanDistance ...@@ -61,8 +61,8 @@ function _find_edge_candidates!(edge_candidates::Vector{EdgeABCEuclideanDistance
return nbr_candidates return nbr_candidates
end end
function _get_edge_index(edge_candidates::Vector{EdgeABCEuclideanDistanceAutomaton}, nbr_candidates::Int, @everywhere function _get_edge_index(edge_candidates::Vector{EdgeABCEuclideanDistanceAutomaton}, nbr_candidates::Int,
detected_event::Bool, tr_nplus1::Transition) detected_event::Bool, tr_nplus1::Transition)
ind_edge = 0 ind_edge = 0
bool_event = detected_event bool_event = detected_event
for i = 1:nbr_candidates for i = 1:nbr_candidates
...@@ -82,11 +82,11 @@ function _get_edge_index(edge_candidates::Vector{EdgeABCEuclideanDistanceAutomat ...@@ -82,11 +82,11 @@ function _get_edge_index(edge_candidates::Vector{EdgeABCEuclideanDistanceAutomat
return (ind_edge, bool_event) return (ind_edge, bool_event)
end end
function next_state!(A::ABCEuclideanDistanceAutomaton, @everywhere function next_state!(A::ABCEuclideanDistanceAutomaton,
ptr_loc_state::Vector{Symbol}, values_state::Vector{Float64}, ptr_time_state::Vector{Float64}, ptr_loc_state::Vector{Symbol}, values_state::Vector{Float64}, ptr_time_state::Vector{Float64},
xnplus1::Vector{Int}, tnplus1::Float64, tr_nplus1::Transition, xnplus1::Vector{Int}, tnplus1::Float64, tr_nplus1::Transition,
xn::Vector{Int}, p::Vector{Float64}, xn::Vector{Int}, p::Vector{Float64},
edge_candidates::Vector{EdgeABCEuclideanDistanceAutomaton}; verbose::Bool = false) edge_candidates::Vector{EdgeABCEuclideanDistanceAutomaton}; verbose::Bool = false)
# En fait d'apres observation de Cosmos, après qu'on ait lu la transition on devrait stop. # En fait d'apres observation de Cosmos, après qu'on ait lu la transition on devrait stop.
detected_event::Bool = false detected_event::Bool = false
turns = 0 turns = 0
...@@ -302,7 +302,7 @@ function create_abc_euclidean_distance_automaton(m::ContinuousTimeModel, timelin ...@@ -302,7 +302,7 @@ function create_abc_euclidean_distance_automaton(m::ContinuousTimeModel, timelin
@everywhere $(update_state!(:l1, :l2, 1))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64) = @everywhere $(update_state!(:l1, :l2, 1))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64) =
(S_values[$(to_idx(:d))] = sqrt(S_values[$(to_idx(:d))]); (S_values[$(to_idx(:d))] = sqrt(S_values[$(to_idx(:d))]);
:l2) :l2)
@everywhere $(check_constraints(:l1, :l2, 2))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64) = @everywhere $(check_constraints(:l1, :l2, 2))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64) =
(S_values[$(to_idx(:d))] > ϵ^2) (S_values[$(to_idx(:d))] > ϵ^2)
@everywhere $(update_state!(:l1, :l2, 2))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64) = @everywhere $(update_state!(:l1, :l2, 2))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64) =
......
...@@ -20,7 +20,11 @@ function _update_values!(values::Vector{Vector{Int}}, times::Vector{Float64}, tr ...@@ -20,7 +20,11 @@ function _update_values!(values::Vector{Vector{Int}}, times::Vector{Float64}, tr
(transitions[idx] = tr_n) (transitions[idx] = tr_n)
end end
function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::Symbol) ########## Code generation methods
# Generates simulate method for a specific ContinuousTimeModel
function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::Symbol;
run_isabsorbing::Bool = false)
return quote return quote
import MarkovProcesses: simulate import MarkovProcesses: simulate
...@@ -72,7 +76,6 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S ...@@ -72,7 +76,6 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S
$(f!)(vec_x, l_t, l_tr, xn, tn, p_sim) $(f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
tn = l_t[1] tn = l_t[1]
isabsorbing = vec_x == xn isabsorbing = vec_x == xn
#isabsorbing = $(isabsorbing)(p_sim,xn)
if isabsorbing || tn > time_bound if isabsorbing || tn > time_bound
break break
end end
...@@ -80,6 +83,10 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S ...@@ -80,6 +83,10 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S
n += 1 n += 1
copyto!(xn, vec_x) copyto!(xn, vec_x)
MarkovProcesses._update_values!(full_values, times, transitions, xn, tn, l_tr[1], i) MarkovProcesses._update_values!(full_values, times, transitions, xn, tn, l_tr[1], i)
if $(run_isabsorbing)
isabsorbing = isabsorbing || $(isabsorbing)(p_sim,xn)
if isabsorbing break end
end
end end
# If simulation ended before the estimation of states # If simulation ended before the estimation of states
if n < estim_min_states if n < estim_min_states
...@@ -101,7 +108,6 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S ...@@ -101,7 +108,6 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S
$(f!)(vec_x, l_t, l_tr, xn, tn, p_sim) $(f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
tn = l_t[1] tn = l_t[1]
isabsorbing = vec_x == xn isabsorbing = vec_x == xn
#isabsorbing = $(isabsorbing)(p_sim,xn)
if isabsorbing || tn > time_bound if isabsorbing || tn > time_bound
i -= 1 i -= 1
break break
...@@ -109,7 +115,11 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S ...@@ -109,7 +115,11 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S
# n+1 values are now n values # n+1 values are now n values
copyto!(xn, vec_x) copyto!(xn, vec_x)
MarkovProcesses._update_values!(full_values, times, transitions, MarkovProcesses._update_values!(full_values, times, transitions,
xn, tn, l_tr[1], estim_min_states+size_tmp+i) xn, tn, l_tr[1], estim_min_states+size_tmp+i)
if $(run_isabsorbing)
isabsorbing = isabsorbing || $(isabsorbing)(p_sim,xn)
if isabsorbing break end
end
end end
# If simulation ended before the end of buffer # If simulation ended before the end of buffer
if i < buffer_size if i < buffer_size
...@@ -129,8 +139,10 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S ...@@ -129,8 +139,10 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S
end end
end end
# Generates simulate method for a specific SynchronizedModel
function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Symbol, function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Symbol,
edge_type::Symbol, f!::Symbol, isabsorbing::Symbol) edge_type::Symbol, f!::Symbol, isabsorbing::Symbol;
run_isabsorbing::Bool = false)
return quote return quote
import MarkovProcesses: simulate, volatile_simulate import MarkovProcesses: simulate, volatile_simulate
...@@ -201,6 +213,10 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym ...@@ -201,6 +213,10 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym
if isacceptedLHA if isacceptedLHA
break break
end end
if $(run_isabsorbing)
isabsorbing = isabsorbing || $(isabsorbing)(p_sim,xn)
if isabsorbing break end
end
end end
# If simulation ended before the estimation of states # If simulation ended before the estimation of states
if n < estim_min_states if n < estim_min_states
...@@ -242,6 +258,10 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym ...@@ -242,6 +258,10 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym
if isacceptedLHA if isacceptedLHA
break break
end end
if $(run_isabsorbing)
isabsorbing = isabsorbing || $(isabsorbing)(p_sim,xn)
if isabsorbing break end
end
end end
# If simulation ended before the end of buffer # If simulation ended before the end of buffer
if i < buffer_size if i < buffer_size
...@@ -311,11 +331,15 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym ...@@ -311,11 +331,15 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym
next_state!(A, ptr_loc_state, values_state, ptr_time_state, next_state!(A, ptr_loc_state, values_state, ptr_time_state,
vec_x, l_t[1], l_tr[1], xn, p_sim, edge_candidates; verbose = verbose) vec_x, l_t[1], l_tr[1], xn, p_sim, edge_candidates; verbose = verbose)
# n+1 values are now n values # n+1 values are now n values
n += 1
copyto!(xn, vec_x) copyto!(xn, vec_x)
tn = l_t[1] tn = l_t[1]
tr_n = l_tr[1] tr_n = l_tr[1]
isacceptedLHA = isaccepted(ptr_loc_state[1], A) isacceptedLHA = isaccepted(ptr_loc_state[1], A)
n += 1 if $(run_isabsorbing)
isabsorbing = isabsorbing || $(isabsorbing)(p_sim,xn)
if isabsorbing break end
end
end end
if isabsorbing && !isacceptedLHA if isabsorbing && !isacceptedLHA
next_state!(A, ptr_loc_state, values_state, ptr_time_state, next_state!(A, ptr_loc_state, values_state, ptr_time_state,
...@@ -328,6 +352,8 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym ...@@ -328,6 +352,8 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym
end end
end end
########## Generic methods
""" """
`volatile_simulate(sm::SynchronizedModel; p, verbose)` `volatile_simulate(sm::SynchronizedModel; p, verbose)`
......
...@@ -24,7 +24,7 @@ nbr_pa = 404 ...@@ -24,7 +24,7 @@ nbr_pa = 404
r = automaton_abc(pm_sync_ER; nbr_particles = nbr_pa) r = automaton_abc(pm_sync_ER; nbr_particles = nbr_pa)
rmprocs(2) rmprocs(workers())
test = size(r.mat_p_end)[1] == pm_sync_ER.df && test = size(r.mat_p_end)[1] == pm_sync_ER.df &&
size(r.mat_p_end)[2] == nbr_pa && size(r.mat_p_end)[2] == nbr_pa &&
......
using MarkovProcesses
using Distributed
addprocs(2)
module_path = get_module_path()
@everywhere module_path = $module_path
@everywhere push!(LOAD_PATH, "$(module_path)/core")
@everywhere using MarkovProcesses
import LinearAlgebra: dot
# ER model
load_model("ER")
set_param!(ER, [:k1, :k2, :k3], [0.2, 40.0, 1.0])
# Observations
timeline = 0.0:0.1:2.0
observations = [vectorize(simulate(ER), :P, timeline)]
epsilon = 0.1 * sqrt(dot(observations[1], observations[1]))
load_automaton("abc_euclidean_distance_automaton")
aut = create_abc_euclidean_distance_automaton(ER, timeline, observations[1], :P)
sync_ER = ER * aut
pm_sync_ER = ParametricModel(sync_ER, (:k1, Uniform(0.0, 100.0)), (:k2, Uniform(0.0, 100.0)))
path_res = "abc_eucl_aut/"
res_abc = automaton_abc(pm_sync_ER; nbr_particles = 100, tolerance = epsilon)
rmprocs(workers())
#=
samples_abc_post = res_abc.mat_p_end
samples_weights = res_abc.weights
histogram2d(samples_abc_post[1,:], samples_abc_post[2,:], weights = samples_weights, normalize = :density)
savefig(path_res * "/histogram.svg")
=#
return true
...@@ -4,7 +4,8 @@ using Test ...@@ -4,7 +4,8 @@ using Test
@testset "ABC SMC and automaton-ABC tests" begin @testset "ABC SMC and automaton-ABC tests" begin
@test include("automaton_abc/R1.jl") @test include("automaton_abc/R1.jl")
@test include("automaton_abc/distributed_R1.jl") @test include("automaton_abc/distributed_R1.jl")
@test include("automaton_abc/abc_euclidean_distance.jl") @test include("automaton_abc/abc_euclidean_distance_automaton.jl")
@test include("automaton_abc/distributed_abc_euclidean_distance_automaton.jl")
#@test test_distributed_R1 #@test test_distributed_R1
end end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment