diff --git a/algorithms/_utils_abc.jl b/algorithms/_utils_abc.jl
index 9e4803459e4b167aec9d370e8cf922db854a7745..7ba9c84af4e895dc32b216b47930bea6c0ae6d04 100644
--- a/algorithms/_utils_abc.jl
+++ b/algorithms/_utils_abc.jl
@@ -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)]
diff --git a/algorithms/abc_smc.jl b/algorithms/abc_smc.jl
index 0640de6cd660d323644ae9ff63f67f7672506c1d..ee22a2d942f48472c8d70b7bc46c2fd6ba1d447a 100644
--- a/algorithms/abc_smc.jl
+++ b/algorithms/abc_smc.jl
@@ -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
diff --git a/automata/abc_euclidean_distance_automaton.jl b/automata/abc_euclidean_distance_automaton.jl
new file mode 100644
index 0000000000000000000000000000000000000000..d1894b97c385313148abb9602d01734f501edd13
--- /dev/null
+++ b/automata/abc_euclidean_distance_automaton.jl
@@ -0,0 +1,354 @@
+
+@everywhere import FunctionWrappers: FunctionWrapper
+@everywhere const ABCCheckConstraintsFunction = FunctionWrapper{Bool,Tuple{Float64,Vector{Float64},Vector{Int},Vector{Float64},Float64}}
+@everywhere const ABCUpdateStateFunction = FunctionWrapper{Symbol,Tuple{Float64,Vector{Float64},Vector{Int},Vector{Float64},Float64}}
+
+@everywhere struct EdgeABCEuclideanDistanceAutomaton <: Edge
+ transitions::TransitionSet
+ check_constraints::ABCCheckConstraintsFunction
+ update_state!::ABCUpdateStateFunction
+end
+
+# Creation of the automaton types
+@everywhere mutable struct ABCEuclideanDistanceAutomaton <: LHA
+ transitions::Vector{Transition}
+ locations::Vector{Location}
+ Λ::Dict{Location,InvariantPredicateFunction}
+ locations_init::Vector{Location}
+ locations_final::Vector{Location}
+ map_var_automaton_idx::Dict{VariableAutomaton,Int} # nvar keys : str_var => idx in values
+ flow::Dict{Location,Vector{Float64}} # output of length nvar
+ map_edges::Dict{Location, Dict{Location,Vector{EdgeABCEuclideanDistanceAutomaton}}}
+ constants::Dict{Symbol,Float64}
+ map_var_model_idx::Dict{VariableModel,Int} # of dim d (of a model)
+ ϵ::Float64
+end
+
+# A push! method implementend by myself because of preallocation of edge_candidates
+function _push_edge!(edge_candidates::Vector{<:EdgeABCEuclideanDistanceAutomaton},
+ edge::EdgeABCEuclideanDistanceAutomaton, nbr_candidates::Int)
+ if nbr_candidates < length(edge_candidates)
+ edge_candidates[nbr_candidates+1] = edge
+ else
+ push!(edge_candidates, edge)
+ end
+end
+
+function _find_edge_candidates!(edge_candidates::Vector{EdgeABCEuclideanDistanceAutomaton},
+ edges_from_current_loc::Dict{Location,Vector{EdgeABCEuclideanDistanceAutomaton}},
+ Λ::Dict{Location,InvariantPredicateFunction},
+ S_time::Float64, S_values::Vector{Float64},
+ x::Vector{Int}, p::Vector{Float64}, ϵ::Float64,
+ only_asynchronous::Bool)
+ nbr_candidates = 0
+ for target_loc in keys(edges_from_current_loc)
+ if !Λ[target_loc](x) continue end
+ for edge in edges_from_current_loc[target_loc]
+ if edge.check_constraints(S_time, S_values, x, p, ϵ)
+ if edge.transitions == nothing
+ _push_edge!(edge_candidates, edge, nbr_candidates)
+ nbr_candidates += 1
+ return nbr_candidates
+ else
+ if !only_asynchronous
+ _push_edge!(edge_candidates, edge, nbr_candidates)
+ nbr_candidates += 1
+ end
+ end
+ end
+ end
+ end
+ return nbr_candidates
+end
+
+function _get_edge_index(edge_candidates::Vector{EdgeABCEuclideanDistanceAutomaton}, nbr_candidates::Int,
+ detected_event::Bool, tr_nplus1::Transition)
+ ind_edge = 0
+ bool_event = detected_event
+ for i = 1:nbr_candidates
+ edge = edge_candidates[i]
+ # Asynchronous edge detection: we fire it
+ if edge.transitions == nothing
+ return (i, detected_event)
+ end
+ # Synchronous detection
+ if !detected_event && tr_nplus1 != nothing
+ if (edge.transitions[1] == :ALL) || (tr_nplus1 in edge.transitions)
+ ind_edge = i
+ bool_event = true
+ end
+ end
+ end
+ return (ind_edge, bool_event)
+end
+
+function next_state!(A::ABCEuclideanDistanceAutomaton,
+ ptr_loc_state::Vector{Symbol}, values_state::Vector{Float64}, ptr_time_state::Vector{Float64},
+ xnplus1::Vector{Int}, tnplus1::Float64, tr_nplus1::Transition,
+ xn::Vector{Int}, p::Vector{Float64},
+ 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.
+ detected_event::Bool = false
+ turns = 0
+ Λ = getfield(A, :Λ)
+ flow = getfield(A, :flow)
+ map_edges = getfield(A, :map_edges)
+ ϵ = A.ϵ
+ if verbose
+ println("##### Begin next_state!")
+ @show xnplus1, tnplus1, tr_nplus1
+ end
+ # First, we check the asynchronous transitions
+ while true
+ turns += 1
+ #edge_candidates = empty!(edge_candidates)
+ edges_from_current_loc = map_edges[ptr_loc_state[1]]
+ # Save all edges that satisfies transition predicate (asynchronous ones)
+ nbr_candidates = _find_edge_candidates!(edge_candidates, edges_from_current_loc, Λ,
+ ptr_time_state[1], values_state, xn, p, ϵ, true)
+ # Search the one we must chose, here the event is nothing because
+ # we're not processing yet the next event
+ ind_edge, detected_event = _get_edge_index(edge_candidates, nbr_candidates, detected_event, nothing)
+ # Update the state with the chosen one (if it exists)
+ # Should be xn here
+ #first_round = false
+ if ind_edge > 0
+ firing_edge = edge_candidates[ind_edge]
+ ptr_loc_state[1] = firing_edge.update_state!(ptr_time_state[1], values_state, xn, p, ϵ)
+ else
+ if verbose println("No edge fired") end
+ break
+ end
+ if verbose
+ @show turns
+ @show edge_candidates
+ @show ind_edge, detected_event, nbr_candidates
+ @show ptr_loc_state[1]
+ @show ptr_time_state[1]
+ @show values_state
+ if turns == 500
+ @warn "We've reached 500 turns"
+ end
+ end
+ # For debug
+ #=
+ if turns > 100
+ println("Number of turns in next_state! is suspicious")
+ @show first_round, detected_event
+ @show length(edge_candidates)
+ )@show tnplus1, tr_nplus1, xnplus1
+ @show edge_candidates
+ error("Unpredicted behavior automaton")
+ end
+ =#
+ end
+ if verbose
+ println("Time flies with the flow...")
+ end
+ # Now time flies according to the flow
+ for i in eachindex(values_state)
+ coeff_deriv = flow[ptr_loc_state[1]][i]
+ if coeff_deriv > 0
+ values_state[i] += coeff_deriv*(tnplus1 - ptr_time_state[1])
+ end
+ end
+ ptr_time_state[1] = tnplus1
+ if verbose
+ @show ptr_loc_state[1]
+ @show ptr_time_state[1]
+ @show values_state
+ end
+ # Now firing an edge according to the event
+ while true
+ turns += 1
+ edges_from_current_loc = map_edges[ptr_loc_state[1]]
+ # Save all edges that satisfies transition predicate (synchronous ones)
+ nbr_candidates = _find_edge_candidates!(edge_candidates, edges_from_current_loc, Λ,
+ ptr_time_state[1], values_state, xnplus1, p, ϵ, false)
+ # Search the one we must chose
+ ind_edge, detected_event = _get_edge_index(edge_candidates, nbr_candidates, detected_event, tr_nplus1)
+ # Update the state with the chosen one (if it exists)
+ if ind_edge > 0
+ firing_edge = edge_candidates[ind_edge]
+ ptr_loc_state[1] = firing_edge.update_state!(ptr_time_state[1], values_state, xnplus1, p, ϵ)
+ end
+ if ind_edge == 0 || detected_event
+ if verbose
+ if detected_event
+ println("Synchronized with $(tr_nplus1)")
+ @show turns
+ @show edge_candidates
+ @show ind_edge, detected_event, nbr_candidates
+ @show detected_event
+ @show ptr_loc_state[1]
+ @show ptr_time_state[1]
+ @show values_state
+ else
+ println("No edge fired")
+ end
+ end
+ break
+ end
+ if verbose
+ @show turns
+ @show edge_candidates
+ @show ind_edge, detected_event, nbr_candidates
+ @show detected_event
+ @show ptr_loc_state[1]
+ @show ptr_time_state[1]
+ @show values_state
+ if turns == 500
+ @warn "We've reached 500 turns"
+ end
+ end
+ # For debug
+ #=
+ if turns > 100
+ println("Number of turns in next_state! is suspicious")
+ @show detected_event
+ @show length(edge_candidates)
+ @show tnplus1, tr_nplus1, xnplus1
+ @show edge_candidates
+ error("Unpredicted behavior automaton")
+ end
+ =#
+ end
+ if verbose
+ println("##### End next_state!")
+ end
+end
+
+function create_abc_euclidean_distance_automaton(m::ContinuousTimeModel, timeline::AbstractVector{Float64}, observations::AbstractVector{Float64}, sym_obs::VariableModel)
+ # Requirements for the automaton
+ @assert sym_obs in m.g "$(sym_obs) is not observed."
+ @assert length(timeline) == length(observations) "Timeline and observations vectors don't have the same length"
+ nbr_observations = length(observations)
+
+ # Automaton types and functions
+ model_name = Symbol(typeof(m))
+ lha_name = :ABCEuclideanDistanceAutomaton
+ edge_type = :EdgeABCEuclideanDistanceAutomaton
+ check_constraints = Symbol("check_constraints_$(lha_name)")
+ update_state! = Symbol("update_state_$(lha_name)!")
+
+ # Locations
+ locations = [:l0, :l1, :l2]
+
+ ## Invariant predicates
+ @everywhere true_inv_predicate(x::Vector{Int}) = true
+ Λ_F = Dict{Location,InvariantPredicateFunction}(:l0 => getfield(Main, :true_inv_predicate), :l1 => getfield(Main, :true_inv_predicate),
+ :l2 => getfield(Main, :true_inv_predicate))
+
+ ## Init and final loc
+ locations_init = [:l0]
+ locations_final = [:l2]
+
+ map_var_automaton_idx = Dict{VariableAutomaton,Int}(:t => 1, :n => 2,
+ :d => 3, :idx => 4)
+ vector_flow = [1.0, 0.0, 0.0, 0.0]
+ flow = Dict{Location,Vector{Float64}}(:l0 => vector_flow,
+ :l1 => vector_flow,
+ :l2 => vector_flow)
+
+ ## Edges
+ idx_obs_var = m.map_var_idx[sym_obs]
+ to_idx(var::Symbol) = map_var_automaton_idx[var]
+
+ id = MarkovProcesses.newid()
+ function check_constraints(from_loc::Location, to_loc::Location, edge_number::Int)
+ return Symbol("check_constraints_$(edge_type)_$(from_loc)$(to_loc)_$(edge_number)_$(model_name)_$(id)")
+ end
+ function update_state!(from_loc::Location, to_loc::Location, edge_number::Int)
+ return Symbol("update_state_$(edge_type)_$(from_loc)$(to_loc)_$(edge_number)_$(model_name)_$(id)!")
+ end
+
+ ## check_constraints & update_state!
+ meta_funcs = quote
+ # l0 loc
+ # l0 => l1
+ #struct $(edge_name(:l0, :l1, 1)) <: $(edge_type) transitions::Union{Nothing,Vector{Symbol}} end
+ @everywhere $(check_constraints(:l0, :l1, 1))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64) = true
+ @everywhere $(update_state!(:l0, :l1, 1))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64) =
+ (S_values[$(to_idx(:n))] = x[$(idx_obs_var)];
+ S_values[$(to_idx(:d))] = 0.0;
+ S_values[$(to_idx(:idx))] = 1.0;
+ :l1)
+
+ # l1 loc
+ # l1 => l1
+ # Defined below
+ #struct $(edge_name(:l1, :l1, 1)) <: $(edge_type) transitions::Union{Nothing,Vector{Symbol}} end
+ @everywhere $(check_constraints(:l1, :l1, 1))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64) =
+ (tml = $(Tuple(timeline));
+ tml_idx = tml[convert(Int, S_values[$(to_idx(:idx))])];
+ S_values[$(to_idx(:t))] >= tml_idx)
+ @everywhere $(update_state!(:l1, :l1, 1))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64) =
+ (y_obs = $(Tuple(observations));
+ y_obs_idx = y_obs[convert(Int, S_values[$(to_idx(:idx))])];
+ S_values[$(to_idx(:d))] += (S_values[$(to_idx(:n))]-y_obs_idx)^2;
+ S_values[$(to_idx(:idx))] += 1.0;
+ :l1)
+
+ #struct $(edge_name(:l1, :l1, 2)) <: $(edge_type) transitions::Union{Nothing,Vector{Symbol}} end
+ @everywhere $(check_constraints(:l1, :l1, 2))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64) = true
+ @everywhere $(update_state!(:l1, :l1, 2))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64) =
+ (S_values[$(to_idx(:n))] = x[$(idx_obs_var)];
+ :l1)
+
+ # l1 => l2
+ #struct $(edge_name(:l1, :l2, 1)) <: $(edge_type) transitions::Union{Nothing,Vector{Symbol}} end
+ @everywhere $(check_constraints(:l1, :l2, 1))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64) =
+ S_values[$(to_idx(:idx))] >= ($nbr_observations + 1)
+ @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))]);
+ :l2)
+
+ @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)
+ @everywhere $(update_state!(:l1, :l2, 2))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}, ϵ::Float64) =
+ (S_values[$(to_idx(:d))] = Inf;
+ :l2)
+ end
+ eval(meta_funcs)
+
+ @eval begin
+ map_edges = Dict{Location,Dict{Location,Vector{$(edge_type)}}}()
+ for loc in $(locations)
+ map_edges[loc] = Dict{Location,Vector{$(edge_type)}}()
+ end
+
+ ## Edges
+ # l0 loc
+ # l0 => l1
+ edge1 = EdgeABCEuclideanDistanceAutomaton(nothing, $(check_constraints(:l0, :l1, 1)), $(update_state!(:l0, :l1, 1)))
+ map_edges[:l0][:l1] = [edge1]
+
+ # l1 loc
+ # l1 => l1
+ edge1 = EdgeABCEuclideanDistanceAutomaton(nothing, $(check_constraints(:l1, :l1, 1)), $(update_state!(:l1, :l1, 1)))
+ edge2 = EdgeABCEuclideanDistanceAutomaton([:ALL], $(check_constraints(:l1, :l1, 2)), $(update_state!(:l1, :l1, 2)))
+ map_edges[:l1][:l1] = [edge1, edge2]
+
+ # l1 => l2
+ edge1 = EdgeABCEuclideanDistanceAutomaton(nothing, $(check_constraints(:l1, :l2, 1)), $(update_state!(:l1, :l2, 1)))
+ edge2 = EdgeABCEuclideanDistanceAutomaton(nothing, $(check_constraints(:l1, :l2, 2)), $(update_state!(:l1, :l2, 2)))
+ map_edges[:l1][:l2] = [edge1,edge2]
+ end
+
+ ## Constants
+ constants = Dict{Symbol,Float64}(:nbr_obs => nbr_observations)
+ for i = 1:nbr_observations
+ constants[Symbol("tml_$(convert(Float64, i))")] = timeline[i]
+ constants[Symbol("y_$(convert(Float64, i))")] = observations[i]
+ end
+
+ # Updating types and simulation methods
+ @everywhere @eval $(MarkovProcesses.generate_code_synchronized_model_type_def(model_name, lha_name))
+ @everywhere @eval $(MarkovProcesses.generate_code_synchronized_simulation(model_name, lha_name, edge_type, m.f!, m.isabsorbing))
+
+ A = ABCEuclideanDistanceAutomaton(m.transitions, locations, Λ_F, locations_init, locations_final,
+ map_var_automaton_idx, flow, map_edges, constants, m.map_var_idx, Inf)
+
+ return A
+end
+
diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index e9f285cf83b2695d64eee87d2231d4b64665db79..9a17dd2f80af7538a4cd538fc943830815167080 100644
--- a/core/MarkovProcesses.jl
+++ b/core/MarkovProcesses.jl
@@ -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
diff --git a/core/model.jl b/core/model.jl
index ef32b376a21b8addc241cab52d2ac0af2c7124f4..92e556af61878fd042d87f6141e6d4c4672a94c5 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -1,7 +1,4 @@
-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)`
diff --git a/tests/automata/euclidean_distance.jl b/tests/automata/euclidean_distance.jl
index 0365a54fabf0f0b015d763e908cd118fed84f104..3bb1dc3a4deb84f8229d6d1d7769e1be3994d868 100644
--- a/tests/automata/euclidean_distance.jl
+++ b/tests/automata/euclidean_distance.jl
@@ -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
diff --git a/tests/automata/euclidean_distance_single.jl b/tests/automata/euclidean_distance_single.jl
index b6f6891614f7fb933d62e09ff5a345462f627e41..dae30d13a66552c755f6798266f58eb683c67450 100644
--- a/tests/automata/euclidean_distance_single.jl
+++ b/tests/automata/euclidean_distance_single.jl
@@ -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
diff --git a/tests/automaton_abc/abc_euclidean_distance_automaton.jl b/tests/automaton_abc/abc_euclidean_distance_automaton.jl
new file mode 100644
index 0000000000000000000000000000000000000000..01a1c03004cedc77a3b6adff83a8c70c23cefd5a
--- /dev/null
+++ b/tests/automaton_abc/abc_euclidean_distance_automaton.jl
@@ -0,0 +1,31 @@
+
+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
+
diff --git a/tests/run_abc_smc.jl b/tests/run_abc_smc.jl
index 8848094b1460ebd2bdd87d1afc5bc88f447130d7..a11f01d46c53b0206d509b75b40c434f5ae636e8 100644
--- a/tests/run_abc_smc.jl
+++ b/tests/run_abc_smc.jl
@@ -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