Commit 43a3f23f authored by Bentriou Mahmoud's avatar Bentriou Mahmoud
Browse files

Meta code generation for next_state! and simulation of synchronized...

Meta code generation for next_state! and simulation of synchronized trajectories. bench/pkg/euclidean_distance*.jl works but segfault with test/automata/euclidean*
parent 5d886fc4
......@@ -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 = m.map_var_idx[sym_obs]
to_idx(var::Symbol) = map_var_automaton_idx[var]
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
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)
@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 = Edge(nothing, getfield(Main, func_name(:cc, :l0, :l1, 1)), getfield(Main, func_name(:us, :l0, :l1, 1)))
edge1 = getfield(Main, $(Meta.quot(edge_name(:l0, :l1, 1))))(nothing)
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)))
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 = Edge(nothing, getfield(Main, func_name(:cc, :l1, :l2, 1)), getfield(Main, func_name(:us, :l1, :l2, 1)))
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
......@@ -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
function generate_code_next_state(lha_name::Symbol, edge_type::Symbol,
check_constraints::Symbol, update_state!::Symbol)
return quote
# 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
end
function generate_next_state_code()
return quote
function _find_edge_candidates!(edge_candidates::Vector{Edge},
edges_from_current_loc::Dict{Location,Vector{Edge}},
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
......
......@@ -132,28 +132,41 @@ function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::S
end
end
function generate_synchronized_simulation_code()
return quote
import MarkovProcesses: simulate
function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Float64}} = nothing,
verbose::Bool = false)
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
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, volatile_simulate
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)