Commit de9fbd7e authored by Bentriou Mahmoud's avatar Bentriou Mahmoud

Add of a new automaton related to ABC.

It computes the euclidean distance and finish the simulation if the
distance is above the tolerance epsilon.
I ha to change some function signature of types in core/model.jl but it
was minor changes.
I add a new test + change some others related to automata.
A minor rewriting core/model for better readability
parent e4859870
......@@ -8,7 +8,7 @@ function _init_abc_automaton!(mat_p_old::Matrix{Float64}, vec_dist::Vector{Float
draw!(vec_p, pm)
mat_p_old[:,i] = vec_p
if l_obs == nothing
S = volatile_simulate(pm, vec_p)
S = volatile_simulate(pm, vec_p; epsilon = Inf)
vec_dist[i] = S[sym_var_aut]
else
l_sim = [simulate(pm, vec_p) for i = 1:length(l_obs)]
......@@ -101,7 +101,7 @@ function _update_param!(mat_p::Matrix{Float64}, vec_dist::Vector{Float64},
continue
end
if l_obs == nothing
S = volatile_simulate(pm, vec_p_prime)
S = volatile_simulate(pm, vec_p_prime; epsilon = epsilon)
dist_sim = S[sym_var_aut]
else
l_sim = [simulate(pm, vec_p_prime) for i = 1:length(l_obs)]
......
......@@ -125,6 +125,7 @@ function _abc_smc(pm::ParametricModel, nbr_particles::Int, tolerance::Float64, a
l_nbr_sim = zeros(Int, nbr_particles)
while (epsilon > last_epsilon) && (current_time - begin_time <= duration_time) && (nbr_tot_sim <= bound_sim)
t += 1
@show pm.m.automaton.ϵ, epsilon
begin_time_ite = time()
@info "Step $t"
# Set new epsilon
......
This diff is collapsed.
......@@ -10,7 +10,9 @@ import Distributed: @everywhere, @distributed
import Distributions: Product, Uniform, Normal
import Distributions: Distribution, Univariate, Continuous, UnivariateDistribution,
MultivariateDistribution, product_distribution
import Distributions: insupport, pdf
import FunctionWrappers: FunctionWrapper
import Random: rand, rand!
import StaticArrays: SVector, @SVector
## Exports
......
import Random: rand, rand!
import Distributions: insupport, pdf
function _resize_trajectory!(values::Vector{Vector{Int}}, times::Vector{Float64},
transitions::Vector{Transition}, size::Int)
for i = eachindex(values) resize!((values[i]), size) end
......@@ -66,7 +63,7 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S
end
return Trajectory(m, values, times, transitions)
end
# Alloc of vectors where we stock n+1 values
# Alloc of vectors where we store n+1 values
vec_x = zeros(Int, getfield(m, :dim_state))
l_t = Float64[0.0]
l_tr = Transition[nothing]
......@@ -74,10 +71,12 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S
for i = 2:estim_min_states
$(f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
tn = l_t[1]
if tn > time_bound || vec_x == xn
isabsorbing = (vec_x == xn)
isabsorbing = vec_x == xn
#isabsorbing = $(isabsorbing)(p_sim,xn)
if isabsorbing || tn > time_bound
break
end
# n+1 values are now n values
n += 1
copyto!(xn, vec_x)
MarkovProcesses._update_values!(full_values, times, transitions, xn, tn, l_tr[1], i)
......@@ -101,15 +100,13 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S
i += 1
$(f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
tn = l_t[1]
if tn > time_bound
i -= 1
break
end
if vec_x == xn
isabsorbing = true
isabsorbing = vec_x == xn
#isabsorbing = $(isabsorbing)(p_sim,xn)
if isabsorbing || tn > time_bound
i -= 1
break
end
# n+1 values are now n values
copyto!(xn, vec_x)
MarkovProcesses._update_values!(full_values, times, transitions,
xn, tn, l_tr[1], estim_min_states+size_tmp+i)
......@@ -189,19 +186,19 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym
# First we fill the allocated array
for i = 2:estim_min_states
$(f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
if l_t[1] > time_bound || vec_x == xn
tn = l_t[1]
isabsorbing = (vec_x == xn)
tn = l_t[1]
isabsorbing = vec_x == xn
if isabsorbing || tn > time_bound
break
end
n += 1
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)
# n+1 values are now n values
n += 1
copyto!(xn, vec_x)
tn = l_t[1]
tr_n = l_tr[1]
MarkovProcesses._update_values!(full_values, times, transitions, xn, tn, tr_n, i)
isacceptedLHA = isaccepted(ptr_loc_state[1], A)
if isabsorbing || isacceptedLHA
if isacceptedLHA
break
end
end
......@@ -229,24 +226,20 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym
while i < buffer_size
i += 1
$(f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
if l_t[1] > time_bound
tn = l_t[1]
i -= 1
break
end
if vec_x == xn
isabsorbing = true
tn = l_t[1]
isabsorbing = vec_x == xn
if isabsorbing || tn > time_bound
i -= 1
break
end
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)
# n+1 values are now n values
copyto!(xn, vec_x)
tn = l_t[1]
tr_n = l_tr[1]
MarkovProcesses._update_values!(full_values, times, transitions,
xn, tn, tr_n, estim_min_states+size_tmp+i)
isacceptedLHA = isaccepted(ptr_loc_state[1], A)
if isabsorbing || isacceptedLHA
if isacceptedLHA
break
end
end
......@@ -259,7 +252,7 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym
end
values = full_values[getfield(m, :_g_idx)]
if isbounded(m) && !isaccepted(ptr_loc_state[1], A)
# Add last value: the convention is the last transition is nothing,
# Add last value: the convention is that if the last transition is nothing,
# the trajectory is bounded
MarkovProcesses._finish_bounded_trajectory!(values, times, transitions, time_bound)
end
......@@ -272,7 +265,9 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym
return SynchronizedTrajectory(S, product, values, times, transitions)
end
function volatile_simulate(m::$(model_name), A::$(lha_name), p_sim::AbstractVector{Float64}, verbose::Bool)
function volatile_simulate(m::$(model_name), A::$(lha_name), p_sim::AbstractVector{Float64},
epsilon::Float64, verbose::Bool)
if $(Meta.quot(lha_name)) == :ABCEuclideanDistanceAutomaton A.ϵ = epsilon end
x0 = getfield(m, :x0)
t0 = getfield(m, :t0)
time_bound = getfield(m, :time_bound)
......@@ -309,12 +304,13 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym
tn = l_t[1]
break
end
if vec_x == xn
isabsorbing = true
isabsorbing = vec_x == xn
if isabsorbing
break
end
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)
# n+1 values are now n values
copyto!(xn, vec_x)
tn = l_t[1]
tr_n = l_tr[1]
......@@ -340,14 +336,14 @@ in order to improve performance.
It returns the last state of the simulation `S::StateLHA` not a trajectory `σ::SynchronizedTrajectory`.
"""
function volatile_simulate(product::SynchronizedModel;
p::Union{Nothing,AbstractVector{Float64}} = nothing, verbose::Bool = false)
p::Union{Nothing,AbstractVector{Float64}} = nothing, epsilon::Float64 = 0.0, verbose::Bool = false)
m = product.m
A = product.automaton
p_sim = getfield(m, :p)
if p != nothing
p_sim = p
end
S = volatile_simulate(m, A, p_sim, verbose)
S = volatile_simulate(m, A, p_sim, epsilon, verbose)
return S
end
......@@ -363,7 +359,6 @@ function simulate(product::SynchronizedModel;
return σ
end
"""
`simulate(pm::ParametricModel, p_prior::AbstractVector{Float64})
......@@ -387,14 +382,10 @@ It returns `S::StateLHA`, not a trajectory.
function volatile_simulate(pm::ParametricModel, p_prior::AbstractVector{Float64};
epsilon::Union{Nothing,Float64} = nothing)
@assert typeof(pm.m) <: SynchronizedModel
# ABC related automata
if @isdefined(EuclideanDistanceABCAutomaton) && typeof(pm.m.automaton) <: EuclideanDistanceABCAutomaton
nothing
end
full_p = copy(get_proba_model(pm).p)
full_p[pm._param_idx] = p_prior
return volatile_simulate(pm.m; p = full_p)
return volatile_simulate(pm.m; p = full_p, epsilon = epsilon)
end
"""
`distribute_mean_value_lha(sm::SynchronizedModel, sym_var::Symbol, nbr_stim::Int)`
......
......@@ -5,6 +5,7 @@ import Distributions: Uniform
load_automaton("euclidean_distance_automaton")
load_automaton("euclidean_distance_automaton_2")
load_automaton("abc_euclidean_distance_automaton")
load_model("SIR")
load_model("ER")
observe_all!(SIR)
......@@ -46,7 +47,18 @@ for i = 1:nbr_sim
else
test2 = true
end
global test_all = test_all && test && test2
sync_SIR = SIR * create_abc_euclidean_distance_automaton(SIR, tml_obs, y_obs, :I)
σ = simulate(sync_SIR)
test3 = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
if !test3
@show test3, euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
global err = σ
global tml = tml_obs
global y = y_obs
global sync_model = sync_SIR
break
end
global test_all = test_all && test && test2 && test3
end
end
......@@ -82,7 +94,18 @@ for i = 1:nbr_sim
else
test2 = true
end
global test_all = test_all && test && test2
sync_ER = ER * create_abc_euclidean_distance_automaton(ER, tml_obs, y_obs, :P)
σ = simulate(sync_ER)
test3 = euclidean_distance(σ, :P, tml_obs, y_obs) == σ.state_lha_end[:d]
if !test3
@show test3, euclidean_distance(σ, :P, tml_obs, y_obs), σ.state_lha_end[:d]
global err = σ
global tml = tml_obs
global y = y_obs
global sync_model = sync_ER
break
end
global test_all = test_all && test && test2 && test3
end
end
......
......@@ -16,12 +16,21 @@ sync_SIR = SIR * aut1
σ = simulate(sync_SIR)
test = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
if !test
@show euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
@show σ
end
if MAKE_SECOND_AUTOMATON_TESTS
load_automaton("euclidean_distance_automaton_2")
aut2 = create_euclidean_distance_automaton_2(SIR, tml_obs, y_obs, :I)
sync_SIR = SIR * aut2
σ = simulate(sync_SIR)
test2 = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
if !test2
@show euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
@show σ
end
else
test2 = true
end
......
using Plots
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)
#=
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,6 +4,7 @@ using Test
@testset "ABC SMC and automaton-ABC tests" begin
@test include("automaton_abc/R1.jl")
@test include("automaton_abc/distributed_R1.jl")
@test include("automaton_abc/abc_euclidean_distance.jl")
#@test test_distributed_R1
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