diff --git a/automata/euclidean_distance_automaton.jl b/automata/euclidean_distance_automaton.jl
index 22d633a427c5feca938599e4708b657ae151b30d..f4c09aaa033e13ee63dba034dd136d57ec257255 100644
--- a/automata/euclidean_distance_automaton.jl
+++ b/automata/euclidean_distance_automaton.jl
@@ -5,13 +5,21 @@ function create_euclidean_distance_automaton(m::ContinuousTimeModel, timeline::A
@assert length(timeline) == length(observations) "Timeline and observations vectors don't have the same length"
nbr_observations = length(observations)
+ # Creation of the automaton types
+ lha_name = :EuclideanDistanceAutomaton
+ edge_type = :EdgeEuclideanDistanceAutomaton
+ check_constraints = Symbol("check_constraints_$(lha_name)")
+ update_state! = Symbol("update_state_$(lha_name)!")
+ @everywhere @eval abstract type $(edge_type) <: Edge end
+ @everywhere @eval $(MarkovProcesses.generate_code_lha_type_def(lha_name, edge_type))
+
# Locations
locations = [:l0, :l1, :l2]
## Invariant predicates
- true_inv_predicate(x::Vector{Int}) = true
- Λ_F = Dict(:l0 => true_inv_predicate, :l1 => true_inv_predicate,
- :l2 => true_inv_predicate)
+ @everywhere true_inv_predicate(x::Vector{Int}) = true
+ Λ_F = Dict{Symbol,Function}(: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]
@@ -25,90 +33,81 @@ function create_euclidean_distance_automaton(m::ContinuousTimeModel, timeline::A
:l2 => vector_flow)
## Edges
- map_edges = Dict{Location, Dict{Location, Vector{Edge}}}()
- for loc in locations
- map_edges[loc] = Dict{Location, Vector{Edge}}()
- end
-
- idx_obs_var = getfield(m, :map_var_idx)[sym_obs]
- idx_var_t = map_var_automaton_idx[:t]
- idx_var_n = map_var_automaton_idx[:n]
- idx_var_d = map_var_automaton_idx[:d]
- idx_var_idx = map_var_automaton_idx[:idx]
-
- nbr_rand = rand(1:1000)
- basename_func = "$(replace(m.name, ' '=>'_'))_$(nbr_rand)"
- basename_func = replace(basename_func, '-'=>'_')
- func_name(type_func::Symbol, from_loc::Location, to_loc::Location, edge_number::Int) =
- Symbol("$(type_func)_eucl_dist_aut_$(basename_func)_$(from_loc)$(to_loc)_$(edge_number)$(type_func == :us ? "!" : "")")
- meta_elementary_functions = quote
+ idx_obs_var = m.map_var_idx[sym_obs]
+ to_idx(var::Symbol) = map_var_automaton_idx[var]
+
+ id = MarkovProcesses.newid()
+ model_name = Symbol(typeof(m))
+ basename_func = "$(model_name)_$(id)"
+ edge_name(from_loc::Location, to_loc::Location, edge_number::Int) =
+ Symbol("Edge_$(lha_name)_$(basename_func)_$(from_loc)$(to_loc)_$(edge_number)")
+
+ ## check_constraints & update_state!
+ meta_elementary_func = quote
# l0 loc
# l0 => l1
- @everywhere $(func_name(:cc, :l0, :l1, 1))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) = true
- @everywhere $(func_name(:us, :l0, :l1, 1))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
- (setindex!(S_values, x[$(idx_obs_var)], $(idx_var_n));
- setindex!(S_values, 0.0, $(idx_var_d));
- setindex!(S_values, 1.0, $(idx_var_idx));
+ struct $(edge_name(:l0, :l1, 1)) <: $(edge_type) transitions::Union{Nothing,Vector{Symbol}} end
+ $(check_constraints)(edge::$(edge_name(:l0, :l1, 1)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) = true
+ $(update_state!)(edge::$(edge_name(:l0, :l1, 1)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
+ (setindex!(S_values, x[$(idx_obs_var)], $(to_idx(:n)));
+ setindex!(S_values, 0.0, $(to_idx(:d)));
+ setindex!(S_values, 1.0, $(to_idx(:idx)));
:l1)
# l1 loc
# l1 => l1
# Defined below
- @everywhere $(func_name(:cc, :l1, :l1, 1))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
+ struct $(edge_name(:l1, :l1, 1)) <: $(edge_type) transitions::Union{Nothing,Vector{Symbol}} end
+ $(check_constraints)(edge::$(edge_name(:l1, :l1, 1)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
(tml = $(Tuple(timeline));
- tml_idx = tml[convert(Int, S_values[$(idx_var_idx)])];
- S_values[$(idx_var_t)] >= tml_idx)
- @everywhere $(func_name(:us, :l1, :l1, 1))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
+ tml_idx = tml[convert(Int, S_values[$(to_idx(:idx))])];
+ S_values[$(to_idx(:t))] >= tml_idx)
+ $(update_state!)(edge::$(edge_name(:l1, :l1, 1)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
(y_obs = $(Tuple(observations));
- y_obs_idx = y_obs[convert(Int, S_values[$(idx_var_idx)])];
- setindex!(S_values, S_values[$(idx_var_d)] + (S_values[$(idx_var_n)] - y_obs_idx)^2,
- $(idx_var_d));
- setindex!(S_values, S_values[$(idx_var_idx)] + 1.0, $(idx_var_idx));
+ y_obs_idx = y_obs[convert(Int, S_values[$(to_idx(:idx))])];
+ setindex!(S_values, S_values[$(to_idx(:d))] + (S_values[$(to_idx(:n))] - y_obs_idx)^2,
+ $(to_idx(:d)));
+ setindex!(S_values, S_values[$(to_idx(:idx))] + 1.0, $(to_idx(:idx)));
:l1)
- @everywhere $(func_name(:cc, :l1, :l1, 2))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) = true
- @everywhere $(func_name(:us, :l1, :l1, 2))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
- (setindex!(S_values, x[$(idx_obs_var)], $(idx_var_n));
+ struct $(edge_name(:l1, :l1, 2)) <: $(edge_type) transitions::Union{Nothing,Vector{Symbol}} end
+ $(check_constraints)(edge::$(edge_name(:l1, :l1, 2)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) = true
+ $(update_state!)(edge::$(edge_name(:l1, :l1, 2)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
+ (setindex!(S_values, x[$(idx_obs_var)], $(to_idx(:n)));
:l1)
# l1 => l2
- @everywhere $(func_name(:cc, :l1, :l2, 1))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
- S_values[$(idx_var_idx)] >= ($nbr_observations + 1)
- @everywhere $(func_name(:us, :l1, :l2, 1))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
- (setindex!(S_values, sqrt(S_values[$(idx_var_d)]), $(idx_var_d));
+ struct $(edge_name(:l1, :l2, 1)) <: $(edge_type) transitions::Union{Nothing,Vector{Symbol}} end
+ $(check_constraints)(edge::$(edge_name(:l1, :l2, 1)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
+ S_values[$(to_idx(:idx))] >= ($nbr_observations + 1)
+ $(update_state!)(edge::$(edge_name(:l1, :l2, 1)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
+ (setindex!(S_values, sqrt(S_values[$(to_idx(:d))]), $(to_idx(:d)));
:l2)
end
- eval(meta_elementary_functions)
-
- # l0 loc
- # l0 => l1
- edge1 = Edge(nothing, getfield(Main, func_name(:cc, :l0, :l1, 1)), getfield(Main, func_name(:us, :l0, :l1, 1)))
- map_edges[:l0][:l1] = [edge1]
-
- # l1 loc
- # l1 => l1
- #=
- edge1 = Edge([:ALL], getfield(Main, func_name(:cc, :l1, :l1, 1)), getfield(Main, func_name(:us, :l1, :l2, 1)))
- map_edges[:l1][:l1] = [edge1]
- for i = 1:nbr_observations
- meta_edge_i = quote
- @everywhere $(func_name(:cc, :l1, :l1, 1+i))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
- S[:t] >= $(timeline[i])
- @everywhere $(func_name(:us, :l1, :l1, 1+i))(S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
- (setindex!(S_values, S[:d] + (S[:n] - $(observations[i]))^2, $(idx_var_d));
- setindex!(S_values, S[:idx] + 1.0, $(idx_var_idx)))
- end
- eval(meta_edge_i)
- push!(map_edges[:l1][:l1], Edge(nothing, getfield(Main, func_name(:cc, :l1, :l1, 1+i)), getfield(Main, func_name(:us, :l1, :l1, 1+i))))
- end
- =#
- edge1 = Edge(nothing, getfield(Main, func_name(:cc, :l1, :l1, 1)), getfield(Main, func_name(:us, :l1, :l1, 1)))
- edge2 = Edge([:ALL], getfield(Main, func_name(:cc, :l1, :l1, 2)), getfield(Main, func_name(:us, :l1, :l1, 2)))
- map_edges[:l1][:l1] = [edge1, edge2]
+ @everywhere eval($(meta_elementary_func))
+
+ @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 = getfield(Main, $(Meta.quot(edge_name(:l0, :l1, 1))))(nothing)
+ map_edges[:l0][:l1] = [edge1]
- # l1 => l2
- edge1 = Edge(nothing, getfield(Main, func_name(:cc, :l1, :l2, 1)), getfield(Main, func_name(:us, :l1, :l2, 1)))
- map_edges[:l1][:l2] = [edge1]
+ # l1 loc
+ # l1 => l1
+ edge1 = getfield(Main, $(Meta.quot(edge_name(:l1, :l1, 1))))(nothing)
+ edge2 = getfield(Main, $(Meta.quot(edge_name(:l1, :l1, 2))))([:ALL])
+ map_edges[:l1][:l1] = [edge1, edge2]
+
+ # l1 => l2
+ edge1 = getfield(Main, $(Meta.quot(edge_name(:l1, :l2, 1))))(nothing)
+ map_edges[:l1][:l2] = [edge1]
+ end
## Constants
constants = Dict{Symbol,Float64}(:nbr_obs => nbr_observations)
@@ -117,10 +116,15 @@ function create_euclidean_distance_automaton(m::ContinuousTimeModel, timeline::A
constants[Symbol("y_$(convert(Float64, i))")] = observations[i]
end
- A = LHA("Euclidean distance", m.transitions, locations, Λ_F, locations_init, locations_final,
- map_var_automaton_idx, flow, map_edges, constants, m.map_var_idx)
+ # Updating next_state!
+ @everywhere @eval $(MarkovProcesses.generate_code_next_state(lha_name, edge_type, check_constraints, update_state!))
+ @everywhere @eval $(MarkovProcesses.generate_code_synchronized_simulation(lha_name, edge_type, m.f!, m.isabsorbing))
+
+ @eval begin
+ A = $(lha_name)($(m.transitions), $(locations), $(Λ_F), $(locations_init), $(locations_final),
+ $(map_var_automaton_idx), $(flow), $(map_edges), $(constants), $(m.map_var_idx))
+ end
+
return A
end
-export create_euclidean_distance_automaton
-
diff --git a/core/lha.jl b/core/lha.jl
index 54f86ca0f6b35b6256b9df4630650c1ef4ac412a..acf19da07fe833d81062124f372d7706c08e3222 100644
--- a/core/lha.jl
+++ b/core/lha.jl
@@ -3,7 +3,7 @@ length_var(A::LHA) = length(getfield(A, :map_var_automaton_idx))
get_value(S::StateLHA, x::Vector{Int}, var::VariableModel) = x[getfield(getfield(S, :A), :map_var_model_idx)[var]]
get_value(A::LHA, x::Vector{Int}, var::VariableModel) = x[getfield(A, :map_var_model_idx)[var]]
-copy(S::StateLHA) = StateLHA(getfield(S, :A), getfield(S, :loc), getfield(S, :values), getfield(S, :time))
+copy(S::StateLHA) = StateLHA(getfield(S, :A), getfield(S, :loc), copy(getfield(S, :values)), getfield(S, :time))
# Not overring getproperty, setproperty to avoid a conversion Symbol => String for the dict key
# From the variable automaton var symbol this function get the index in S.values
@@ -18,7 +18,7 @@ setindex!(S::StateLHA, val::Int, var::VariableAutomaton) = S[var] = convert(Floa
setindex!(S::StateLHA, val::Bool, var::VariableAutomaton) = S[var] = convert(Float64, val)
function Base.show(io::IO, A::LHA)
- print(io, "$(A.name) automaton (LHA)\n")
+ print(io, "$(Symbol(typeof(A))) automaton (LHA)\n")
print(io, "- initial locations : $(join(A.locations_init,','))\n")
print(io, "- final locations : $(join(A.locations_final,','))\n")
print(io, "- labeling prop Λ :\n")
@@ -60,10 +60,6 @@ end
function Base.show(io::IO, E::Edge)
print(io, "($(typeof(E)): ")
print(io, (E.transitions == nothing) ? "Asynchronous #)" : ("Synchronized with " * join(E.transitions,',') * ")"))
- #=print(io, Symbol(E.check_constraints))
- print(io, ", ")
- print(io, Symbol(E.update_state!))
- print(io, ")")=#
end
function Base.copyto!(Sdest::StateLHA, Ssrc::StateLHA)
@@ -79,6 +75,7 @@ end
# of one of the LHA fields
isaccepted(S::StateLHA) = (getfield(S, :loc) in getfield(getfield(S, :A), :locations_final))
+isaccepted(loc::Symbol, A::LHA) = loc in A.locations_final
# Methods for synchronize / read the trajectory
function init_state(A::LHA, x0::Vector{Int}, t0::Float64)
@@ -92,19 +89,21 @@ function init_state(A::LHA, x0::Vector{Int}, t0::Float64)
return S0
end
-# A push! method implementend by myself because of preallocation of edge_candidates
-function _push_edge!(edge_candidates::Vector{Edge}, edge::Edge, nbr_candidates::Int)
- if nbr_candidates < 2
- @inbounds edge_candidates[nbr_candidates+1] = edge
- else
- push!(edge_candidates, edge)
- end
-end
-
-function generate_next_state_code()
+function generate_code_next_state(lha_name::Symbol, edge_type::Symbol,
+ check_constraints::Symbol, update_state!::Symbol)
+
return quote
- function _find_edge_candidates!(edge_candidates::Vector{Edge},
- edges_from_current_loc::Dict{Location,Vector{Edge}},
+ # A push! method implementend by myself because of preallocation of edge_candidates
+ function _push_edge!(edge_candidates::Vector{<:$(edge_type)}, edge::$(edge_type), nbr_candidates::Int)
+ if nbr_candidates < length(edge_candidates)
+ @inbounds edge_candidates[nbr_candidates+1] = edge
+ else
+ push!(edge_candidates, edge)
+ end
+ end
+
+ function _find_edge_candidates!(edge_candidates::Vector{$(edge_type)},
+ edges_from_current_loc::Dict{Location,Vector{$(edge_type)}},
Λ::Dict{Location,Function},
S_time::Float64, S_values::Vector{Float64},
x::Vector{Int}, p::Vector{Float64},
@@ -113,14 +112,14 @@ function generate_next_state_code()
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 check_constraints(edge, S_time, S_values, x, p)
+ if $(check_constraints)(edge, S_time, S_values, x, p)
if getfield(edge, :transitions) == nothing
- MarkovProcesses._push_edge!(edge_candidates, edge, nbr_candidates)
+ _push_edge!(edge_candidates, edge, nbr_candidates)
nbr_candidates += 1
return nbr_candidates
else
if !only_asynchronous
- MarkovProcesses._push_edge!(edge_candidates, edge, nbr_candidates)
+ _push_edge!(edge_candidates, edge, nbr_candidates)
nbr_candidates += 1
end
end
@@ -130,7 +129,7 @@ function generate_next_state_code()
return nbr_candidates
end
- function _get_edge_index(edge_candidates::Vector{Edge}, nbr_candidates::Int,
+ function _get_edge_index(edge_candidates::Vector{$(edge_type)}, nbr_candidates::Int,
detected_event::Bool, tr_nplus1::Transition)
ind_edge = 0
bool_event = detected_event
@@ -152,16 +151,14 @@ function generate_next_state_code()
return (ind_edge, bool_event)
end
- function next_state!(Snplus1::StateLHA, A::LHA,
+ function next_state!(A::$(lha_name),
+ ptr_loc_state::Vector{Symbol}, values_state::Vector{Float64}, ptr_time_state::Vector{Float64},
xnplus1::Vector{Int}, tnplus1::Float64, tr_nplus1::Transition,
- Sn::StateLHA, xn::Vector{Int}, p::Vector{Float64},
- edge_candidates::Vector{Edge}; verbose::Bool = false)
+ xn::Vector{Int}, p::Vector{Float64},
+ edge_candidates::Vector{$(edge_type)}; 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
- current_values = getfield(Snplus1, :values)
- current_time = getfield(Snplus1, :time)
- current_loc = getfield(Snplus1, :loc)
Λ = getfield(A, :Λ)
flow = getfield(A, :flow)
map_edges = getfield(A, :map_edges)
@@ -176,10 +173,10 @@ function generate_next_state_code()
while true
turns += 1
#edge_candidates = empty!(edge_candidates)
- edges_from_current_loc = map_edges[current_loc]
+ 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, Λ,
- current_time, current_values, xn, p, true)
+ 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)
@@ -188,7 +185,7 @@ function generate_next_state_code()
#first_round = false
if ind_edge > 0
firing_edge = edge_candidates[ind_edge]
- current_loc = getfield(Main, :update_state!)(firing_edge, current_time, current_values, xn, p)
+ ptr_loc_state[1] = $(update_state!)(firing_edge, ptr_time_state[1], values_state, xn, p)
else
if verbose println("No edge fired") end
break
@@ -197,9 +194,9 @@ function generate_next_state_code()
@show turns
@show edge_candidates
@show ind_edge, detected_event, nbr_candidates
- @show current_loc
- @show current_time
- @show current_values
+ @show ptr_loc_state[1]
+ @show ptr_time_state[1]
+ @show values_state
if turns == 500
@warn "We've reached 500 turns"
end
@@ -212,9 +209,6 @@ function generate_next_state_code()
@show length(edge_candidates)
@show tnplus1, tr_nplus1, xnplus1
@show edge_candidates
- for edge in edge_candidates
- @show getfield(edge, :check_constraints)(time_S, values_S, xn, p)
- end
error("Unpredicted behavior automaton")
end
=#
@@ -223,31 +217,31 @@ function generate_next_state_code()
println("Time flies with the flow...")
end
# Now time flies according to the flow
- for i in eachindex(current_values)
- @inbounds coeff_deriv = flow[current_loc][i]
+ for i in eachindex(values_state)
+ @inbounds coeff_deriv = flow[ptr_loc_state[1]][i]
if coeff_deriv > 0
- @inbounds current_values[i] += coeff_deriv*(tnplus1 - current_time)
+ @inbounds values_state[i] += coeff_deriv*(tnplus1 - ptr_time_state[1])
end
end
- current_time = tnplus1
+ ptr_time_state[1] = tnplus1
if verbose
- @show current_loc
- @show current_time
- @show current_values
+ @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[current_loc]
+ 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, Λ,
- current_time, current_values, xnplus1, p, false)
+ 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]
- current_loc = getfield(Main, :update_state!)(firing_edge, current_time, current_values, xnplus1, p)
+ ptr_loc_state[1] = $(update_state!)(firing_edge, ptr_time_state[1], values_state, xnplus1, p)
end
if ind_edge == 0 || detected_event
if verbose
@@ -257,9 +251,9 @@ function generate_next_state_code()
@show edge_candidates
@show ind_edge, detected_event, nbr_candidates
@show detected_event
- @show current_loc
- @show current_time
- @show current_values
+ @show ptr_loc_state[1]
+ @show ptr_time_state[1]
+ @show values_state
else
println("No edge fired")
end
@@ -271,9 +265,9 @@ function generate_next_state_code()
@show edge_candidates
@show ind_edge, detected_event, nbr_candidates
@show detected_event
- @show current_loc
- @show current_time
- @show current_values
+ @show ptr_loc_state[1]
+ @show ptr_time_state[1]
+ @show values_state
if turns == 500
@warn "We've reached 500 turns"
end
@@ -286,17 +280,15 @@ function generate_next_state_code()
@show length(edge_candidates)
@show tnplus1, tr_nplus1, xnplus1
@show edge_candidates
- for edge in edge_candidates
- @show getfield(edge, :check_constraints)(time_S, values_S, x, p)
- end
error("Unpredicted behavior automaton")
end
=#
end
- setfield!(Snplus1, :loc, current_loc)
- setfield!(Snplus1, :time, current_time)
+ #=
+ setfield!(Snplus1, :loc, ptr_loc_state[1])
+ setfield!(Snplus1, :time, ptr_time_state[1])
+ =#
if verbose
- @show Snplus1
println("##### End next_state!")
end
end
diff --git a/core/model.jl b/core/model.jl
index 5d2b0220bf670799207f611c50960d4811ff5399..e56f9d50d764698ea72df844eb75ff03ff681e0e 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -132,28 +132,41 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S
end
end
-function generate_synchronized_simulation_code()
+function simulate(product::SynchronizedModel;
+ p::Union{Nothing,AbstractVector{Float64}} = nothing, verbose::Bool = false)
+ m = getfield(product, :m)
+ A = getfield(product, :automaton)
+ p_sim = getfield(m, :p)
+ if p != nothing
+ p_sim = p
+ end
+ return simulate(m, A, product, p_sim, verbose)
+end
+
+function volatile_simulate(product::SynchronizedModel;
+ p::Union{Nothing,AbstractVector{Float64}} = nothing, verbose::Bool = false)
+ m = product.m
+ A = product.automaton
+ p_sim = getfield(m, :p)
+ if p != nothing
+ p_sim = p
+ end
+ return volatile_simulate(m, A, p_sim, verbose)
+end
+
+function generate_code_synchronized_simulation(lha_name::Symbol, edge_type::Symbol, f!::Symbol, isabsorbing::Symbol)
return quote
- import MarkovProcesses: simulate
+ import MarkovProcesses: simulate, volatile_simulate
- function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Float64}} = nothing,
- verbose::Bool = false)
- m = getfield(product, :m)
- A = getfield(product, :automaton)
- p_sim = getfield(m, :p)
- func_f! = getfield(m, :f!)
- func_isabsorbing = getfield(m, :isabsorbing)
- if p != nothing
- p_sim = p
- end
+ function simulate(m::ContinuousTimeModel, A::$(lha_name), product::SynchronizedModel,
+ p_sim::AbstractVector{Float64}, verbose::Bool)
x0 = getfield(m, :x0)
t0 = getfield(m, :t0)
time_bound = getfield(m, :time_bound)
- S0 = init_state(A, x0, t0)
buffer_size = getfield(m, :buffer_size)
estim_min_states = getfield(m, :estim_min_states)
- edge_candidates = Vector{Edge}(undef, 2)
+ edge_candidates = Vector{$(edge_type)}(undef, 2)
# First alloc
full_values = Vector{Vector{Int}}(undef, getfield(m, :dim_state))
for i = eachindex(full_values) full_values[i] = zeros(Int, estim_min_states) end
@@ -163,19 +176,22 @@ function generate_synchronized_simulation_code()
for i = eachindex(full_values) full_values[i][1] = x0[i] end
times[1] = t0
transitions[1] = nothing
+ S = init_state(A, x0, t0)
+ ptr_loc_state = [S.loc]
+ values_state = S.values
+ ptr_time_state = [S.time]
# Values at time n
n = 1
xn = copy(x0)
tn = copy(t0)
- Sn = copy(S0)
- next_state!(Sn, A, x0, t0, nothing, S0, x0, p_sim, edge_candidates; verbose = verbose)
- isabsorbing::Bool = func_isabsorbing(p_sim,xn)
- isacceptedLHA::Bool = isaccepted(Sn)
+ next_state!(A, ptr_loc_state, values_state, ptr_time_state,
+ x0, t0, nothing, x0, p_sim, edge_candidates; verbose = verbose)
+ isabsorbing::Bool = $(isabsorbing)(p_sim,xn)
+ isacceptedLHA::Bool = isaccepted(ptr_loc_state[1], A)
# Alloc of vectors where we stock n+1 values
vec_x = zeros(Int, getfield(m, :dim_state))
l_t = Float64[0.0]
l_tr = Transition[nothing]
- Snplus1 = deepcopy(Sn)
# If x0 is absorbing
if isabsorbing || isacceptedLHA
MarkovProcesses._resize_trajectory!(full_values, times, transitions, 1)
@@ -184,27 +200,28 @@ function generate_synchronized_simulation_code()
MarkovProcesses._finish_bounded_trajectory!(values, times, transitions, time_bound)
end
if isabsorbing && !isacceptedLHA
- next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
- copyto!(Sn, Snplus1)
+ next_state!(A, ptr_loc_state, values_state, ptr_time_state,
+ xn, time_bound, nothing, xn, p_sim, edge_candidates; verbose = verbose)
end
- return SynchronizedTrajectory(Sn, product, values, times, transitions)
+ setfield!(S, :loc, ptr_loc_state[1])
+ setfield!(S, :time, ptr_time_state[1])
+ return SynchronizedTrajectory(S, product, values, times, transitions)
end
# First we fill the allocated array
for i = 2:estim_min_states
- func_f!(vec_x, l_t, l_tr, xn, tn, p_sim)
+ $(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)
break
end
n += 1
- next_state!(Snplus1, A, vec_x, l_t[1], l_tr[1], Sn, xn, p_sim, edge_candidates; verbose = verbose)
+ 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)
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)
- copyto!(Sn, Snplus1)
- isacceptedLHA = isaccepted(Snplus1)
+ isacceptedLHA = isaccepted(ptr_loc_state[1], A)
if isabsorbing ||Â isacceptedLHA
break
end
@@ -217,10 +234,12 @@ function generate_synchronized_simulation_code()
MarkovProcesses._finish_bounded_trajectory!(values, times, transitions, time_bound)
end
if isabsorbing && !isacceptedLHA
- next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
- copyto!(Sn, Snplus1)
+ next_state!(A, ptr_loc_state, values_state, ptr_time_state,
+ xn, time_bound, nothing, xn, p_sim, edge_candidates; verbose = verbose)
end
- return SynchronizedTrajectory(Sn, product, values, times, transitions)
+ setfield!(S, :loc, ptr_loc_state[1])
+ setfield!(S, :time, ptr_time_state[1])
+ return SynchronizedTrajectory(S, product, values, times, transitions)
end
# Otherwise, buffering system
size_tmp = 0
@@ -230,7 +249,7 @@ function generate_synchronized_simulation_code()
i = 0
while i < buffer_size
i += 1
- func_f!(vec_x, l_t, l_tr, xn, tn, p_sim)
+ $(f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
if l_t[1] > time_bound
tn = l_t[1]
i -= 1
@@ -241,14 +260,13 @@ function generate_synchronized_simulation_code()
i -= 1
break
end
- next_state!(Snplus1, A, vec_x, l_t[1], l_tr[1], Sn, xn, p_sim, edge_candidates; verbose = verbose)
+ 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)
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)
- copyto!(Sn, Snplus1)
- isacceptedLHA = isaccepted(Snplus1)
+ xn, tn, tr_n, estim_min_states+size_tmp+i)
+ isacceptedLHA = isaccepted(ptr_loc_state[1], A)
if isabsorbing || isacceptedLHA
break
end
@@ -261,17 +279,20 @@ function generate_synchronized_simulation_code()
n += i
end
values = full_values[getfield(m, :_g_idx)]
- if isbounded(m) && !isaccepted(Sn)
+ if isbounded(m) && !isaccepted(ptr_loc_state[1], A)
# Add last value: the convention is the last transition is nothing,
# the trajectory is bounded
MarkovProcesses._finish_bounded_trajectory!(values, times, transitions, time_bound)
end
if isabsorbing && !isacceptedLHA
- next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
- copyto!(Sn, Snplus1)
+ next_state!(A, ptr_loc_state, values_state, ptr_time_state,
+ xn, time_bound, nothing, xn, p_sim, edge_candidates; verbose = verbose)
end
- return SynchronizedTrajectory(Sn, product, values, times, transitions)
+ setfield!(S, :loc, ptr_loc_state[1])
+ setfield!(S, :time, ptr_time_state[1])
+ return SynchronizedTrajectory(S, product, values, times, transitions)
end
+
"""
`volatile_simulate(sm::SynchronizedModel; p, verbose)`
@@ -279,44 +300,39 @@ function generate_synchronized_simulation_code()
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)
- m = product.m
- A = product.automaton
- p_sim = getfield(m, :p)
- func_f! = getfield(m, :f!)
- func_isabsorbing = getfield(m, :isabsorbing)
- if p != nothing
- p_sim = p
- end
+ function volatile_simulate(m::ContinuousTimeModel, A::$(lha_name), p_sim::AbstractVector{Float64}, verbose::Bool)
x0 = getfield(m, :x0)
t0 = getfield(m, :t0)
time_bound = getfield(m, :time_bound)
- S0 = init_state(A, x0, t0)
- edge_candidates = Vector{Edge}(undef, 2)
+ S = init_state(A, x0, t0)
+ ptr_loc_state = [S.loc]
+ values_state = S.values
+ ptr_time_state = [S.time]
+ edge_candidates = Vector{$(edge_type)}(undef, 2)
# Values at time n
n = 1
xn = copy(x0)
tn = copy(t0)
- Sn = copy(S0)
- next_state!(Sn, A, x0, t0, nothing, S0, x0, p_sim, edge_candidates; verbose = verbose)
- isabsorbing::Bool = func_isabsorbing(p_sim,xn)
- isacceptedLHA::Bool = isaccepted(Sn)
+ next_state!(A, ptr_loc_state, values_state, ptr_time_state,
+ x0, t0, nothing, x0, p_sim, edge_candidates; verbose = verbose)
+ isabsorbing::Bool = $(isabsorbing)(p_sim,xn)
+ isacceptedLHA::Bool = isaccepted(ptr_loc_state[1], A)
# Alloc of vectors where we stock n+1 values
vec_x = zeros(Int, getfield(m, :dim_state))
l_t = Float64[0.0]
l_tr = Transition[nothing]
- Snplus1 = deepcopy(Sn)
# If x0 is absorbing
if isabsorbing || isacceptedLHA
if !isacceptedLHA
- next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
- copyto!(Sn, Snplus1)
+ next_state!(A, ptr_loc_state, values_state, ptr_time_state,
+ xn, time_bound, nothing, xn, p_sim, edge_candidates; verbose = verbose)
end
- return Sn
+ setfield!(S, :loc, ptr_loc_state[1])
+ setfield!(S, :time, ptr_time_state[1])
+ return S
end
while !isabsorbing &&Â tn <= time_bound && !isacceptedLHA
- func_f!(vec_x, l_t, l_tr, xn, tn, p_sim)
+ $(f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
if l_t[1] > time_bound
tn = l_t[1]
break
@@ -325,19 +341,21 @@ function generate_synchronized_simulation_code()
isabsorbing = true
break
end
- next_state!(Snplus1, A, vec_x, l_t[1], l_tr[1], Sn, xn, p_sim, edge_candidates; verbose = verbose)
+ 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)
copyto!(xn, vec_x)
tn = l_t[1]
tr_n = l_tr[1]
- copyto!(Sn, Snplus1)
- isacceptedLHA = isaccepted(Snplus1)
+ isacceptedLHA = isaccepted(ptr_loc_state[1], A)
n += 1
end
if isabsorbing && !isacceptedLHA
- next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
- copyto!(Sn, Snplus1)
+ next_state!(A, ptr_loc_state, values_state, ptr_time_state,
+ xn, time_bound, nothing, xn, p_sim, edge_candidates; verbose = verbose)
end
- return Sn
+ setfield!(S, :loc, ptr_loc_state[1])
+ setfield!(S, :time, ptr_time_state[1])
+ return S
end
end
end
@@ -471,7 +489,7 @@ function check_consistency(m::ContinuousTimeModel)
@assert length(m.g) <= m.dim_state
@assert length(m._g_idx) == length(m.g)
@assert m.buffer_size >= 0
- @assert typeof(getfield(Main, m.isabsorbing)(m.p, m.x0)) == Bool
+ #@assert typeof(getfield(Main, m.isabsorbing)(m.p, m.x0)) == Bool
return true
end