From 21c72d7eb86b0ef47ff45529eea976f4d8cbdfee Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Fri, 20 Nov 2020 21:19:46 +0100
Subject: [PATCH] Add of several useful functions for data access for
Trajectory object.
First implementation of LHA.
Synchronized simulation not yet tested.
Reading of trajectories with LHA seems to work with one test basted on
simulations that checks that the last value is well read.
---
automata/automaton_F.jl | 95 +++++++++++
core/MarkovProcesses.jl | 12 +-
core/lha.jl | 161 ++++++++++++++++++
core/model.jl | 72 +++++++-
core/trajectory.jl | 51 ++++--
tests/automata/automaton_F.jl | 12 ++
.../automata/read_trajectory_last_state_F.jl | 32 ++++
tests/run_all.jl | 1 +
tests/run_automata.jl | 7 +
tests/run_bench.jl | 16 ++
10 files changed, 436 insertions(+), 23 deletions(-)
create mode 100644 automata/automaton_F.jl
create mode 100644 core/lha.jl
create mode 100644 tests/automata/automaton_F.jl
create mode 100644 tests/automata/read_trajectory_last_state_F.jl
create mode 100644 tests/run_automata.jl
create mode 100644 tests/run_bench.jl
diff --git a/automata/automaton_F.jl b/automata/automaton_F.jl
new file mode 100644
index 0000000..38c2c3e
--- /dev/null
+++ b/automata/automaton_F.jl
@@ -0,0 +1,95 @@
+
+function create_automaton_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1::Float64, t2::Float64, str_obs::String)
+ # Locations
+ l_loc = ["l0", "l1", "l2", "l3"]
+
+ ## Invariant predicates
+ true_inv_predicate = (A::LHA, S:: StateLHA) -> return true
+ Λ_F = Dict("l0" => true_inv_predicate, "l1" => true_inv_predicate,
+ "l2" => true_inv_predicate, "l3" => true_inv_predicate)
+
+ ## Init and final loc
+ l_loc_init = ["l0"]
+ l_loc_final = ["l2"]
+
+ #S.n <=> S.l_var[A.map_var_automaton_idx["n"]]
+ #P <=> xn[map_var_model_idx[l_ctes[str_O]] with str_O = "P". On stock str_O dans l_ctes
+
+ ## Map of automaton variables
+ map_var_automaton_idx = Dict{VariableAutomaton,Int}("n" => 1, "d" => 2)
+
+ ## Flow of variables
+ l_flow = Dict{VariableAutomaton,Vector{Float64}}("l0" => [0.0,0.0],
+ "l1" => [0.0,0.0],
+ "l2" => [0.0,0.0],
+ "l3" => [0.0,0.0])
+
+ ## Edges
+ map_edges = Dict{Tuple{Location,Location}, Vector{Edge}}()
+
+ # l0 loc : we construct the edges of the form l0 => (..)
+ # "cc" as check_constraints
+ tuple = ("l0", "l1")
+ cc_aut_F_l0l1_1(A::LHA, S::StateLHA) = true
+ us_aut_F_l0l1_1!(A::LHA, S::StateLHA, x::SubArray{Int,1}) =
+ (S.loc = "l1"; S["d"] = Inf; S["n"] = get_value(A, x, str_obs))
+ edge1 = Edge([nothing], cc_aut_F_l0l1_1, us_aut_F_l0l1_1!)
+ map_edges[tuple] = [edge1]
+
+ # l1 loc : we construct the edges of the form l1 => (..)
+ tuple = ("l1", "l2")
+ cc_aut_F_l1l2_1(A::LHA, S::StateLHA) =
+ (A.l_ctes["x1"] <= S["n"] <= A.l_ctes["x2"]) && (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"])
+ us_aut_F_l1l2_1!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = (S["d"] = 0; S.loc = "l2")
+ edge1 = Edge([nothing], cc_aut_F_l1l2_1, us_aut_F_l1l2_1!)
+
+ cc_aut_F_l1l2_2(A::LHA, S::StateLHA) = (S["d"] > 0 && S.time > A.l_ctes["t2"])
+ us_aut_F_l1l2_2!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = (S.loc = "l2")
+ edge2 = Edge([nothing], cc_aut_F_l1l2_2, us_aut_F_l1l2_2!)
+
+ cc_aut_F_l1l2_3(A::LHA, S::StateLHA) = (S["d"] == 0 && S.time >= A.l_ctes["t1"])
+ us_aut_F_l1l2_3!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = (S.loc = "l2")
+ edge3 = Edge([nothing], cc_aut_F_l1l2_3, us_aut_F_l1l2_3!)
+
+ map_edges[tuple] = [edge1, edge2, edge3]
+
+ tuple = ("l1", "l3")
+ cc_aut_F_l1l3_1(A::LHA, S::StateLHA) = (A.l_ctes["x1"] <= S["n"] <= A.l_ctes["x2"])
+ us_aut_F_l1l3_1!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = (S["d"] = 0; S.loc = "l3")
+ edge1 = Edge([nothing], cc_aut_F_l1l3_1, us_aut_F_l1l3_1!)
+
+ cc_aut_F_l1l3_2(A::LHA, S::StateLHA) =
+ (A.l_ctes["x1"] > S["n"] || S["n"] > A.l_ctes["x2"]) && (S.time < A.l_ctes["t1"])
+ us_aut_F_l1l3_2!(A::LHA, S::StateLHA, x::SubArray{Int,1}) =
+ (S["d"] = min(sqrt((S.time - A.l_ctes["t1"])^2 + (S["n"] - A.l_ctes["x2"])^2),
+ sqrt((S.time - A.l_ctes["t1"])^2 + (S["n"] - A.l_ctes["x1"])^2)); S.loc = "l3")
+ edge2 = Edge([nothing], cc_aut_F_l1l3_2, us_aut_F_l1l3_2!)
+
+ cc_aut_F_l1l3_3(A::LHA, S::StateLHA) =
+ (A.l_ctes["x1"] > S["n"] || S["n"] > A.l_ctes["x2"]) && (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"])
+ us_aut_F_l1l3_3!(A::LHA, S::StateLHA, x::SubArray{Int,1}) =
+ (S["d"] = min(S["d"], min(abs(S["n"] - A.l_ctes["x1"]), abs(S["n"] - A.l_ctes["x2"]))); S.loc = "l3")
+ edge3 = Edge([nothing], cc_aut_F_l1l3_3, us_aut_F_l1l3_3!)
+ map_edges[tuple] = [edge1, edge2, edge3]
+
+ # l3 loc
+ tuple = ("l3", "l1")
+ cc_aut_F_l3l1_1(A::LHA, S::StateLHA) = true
+ us_aut_F_l3l1_1!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = (S["n"] = get_value(A, x, str_obs); S.loc = "l1")
+ edge1 = Edge(["ALL"], cc_aut_F_l3l1_1, us_aut_F_l3l1_1!)
+ tuple = ("l3", "l2")
+ cc_aut_F_l3l2_1(A::LHA, S::StateLHA) = (S.time >= A.l_ctes["t2"])
+ us_aut_F_l3l2_1!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = (S["n"] = get_value(A, x, str_obs); S.loc = "l2")
+ edge2 = Edge([nothing], cc_aut_F_l3l2_1, us_aut_F_l3l2_1!)
+ map_edges[tuple] = [edge1, edge2]
+
+ ## Constants
+ l_ctes = Dict{String,Float64}("x1" => x1, "x2" => x2, "t1" => t1, "t2" => t2)
+
+ A = LHA(m.l_transitions, l_loc, Λ_F, l_loc_init, l_loc_final,
+ map_var_automaton_idx, l_flow, map_edges, l_ctes, m.map_var_idx)
+ return A
+end
+
+export create_automaton_F
+
diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index e9d2614..dfb23e6 100644
--- a/core/MarkovProcesses.jl
+++ b/core/MarkovProcesses.jl
@@ -1,6 +1,7 @@
module MarkovProcesses
-import Base: +, -, getfield, getindex, ImmutableDict
+import Base: +, -, *
+import Base: copy, getfield, getindex, lastindex, setindex!, getproperty, setproperty!
export Model, ContinuousTimeModel, DiscreteTimeModel
export simulate, set_param!, get_param, set_observed_var!
@@ -10,8 +11,15 @@ include("model.jl")
export Observations, AbstractTrajectory, Trajectory
export +, -, δ, dist_lp
-export get_obs_var, length_states, length_obs_var, is_bounded, times, transitions
+export get_obs_var, length_states, length_obs_var, get_state_from_time
+export is_bounded, times, transitions
+export check_consistency, is_steadystate
include("trajectory.jl")
+export LHA, StateLHA, Edge
+export init_state, next_state!, read_trajectory
+export load_automaton, get_index, get_value, length_var
+include("lha.jl")
+
end
diff --git a/core/lha.jl b/core/lha.jl
new file mode 100644
index 0000000..c3a1359
--- /dev/null
+++ b/core/lha.jl
@@ -0,0 +1,161 @@
+
+const Location = String
+const VariableAutomaton = String
+
+struct Edge
+ transitions::Vector{Transition}
+ check_constraints::Function
+ update_state!::Function
+end
+
+struct LHA
+ l_transitions::Vector{Transition}
+ l_loc::Vector{Location}
+ Λ::Dict{Location,Function}
+ l_loc_init::Vector{Location}
+ l_loc_final::Vector{Location}
+ map_var_automaton_idx::Dict{VariableAutomaton,Int} # nvar keys : str_var => idx in l_var
+ l_flow::Dict{Location,Vector{Float64}} # output of length nvar
+ map_edges::Dict{Tuple{Location,Location},Vector{Edge}}
+ l_ctes::Dict{String,Float64}
+ map_var_model_idx::Dict{String,Int} # of dim d (of a model)
+end
+
+LHA(A::LHA, map_var::Dict{String,Int}) = LHA(A.l_transitions, A.l_loc, A.Λ,
+ A.l_loc_init, A.l_loc_final, A.map_var_automaton_idx, A.l_flow,
+ A.map_edges, A.l_ctes, map_var)
+length_var(A::LHA) = length(A.map_var_automaton_idx)
+get_value(A::LHA, x::SubArray{Int,1}, var::String) = x[A.map_var_model_idx[var]]
+load_automaton(automaton::String) = include(get_module_path() * "/automata/$(automaton).jl")
+
+mutable struct StateLHA
+ A::LHA
+ loc::Location
+ l_var::Vector{Float64}
+ time::Float64
+end
+
+copy(S::StateLHA) = StateLHA(S.A, S.loc, S.l_var, S.time)
+# Not overring getproperty, setproperty to avoid a conversion Symbol => String for the dict key
+getindex(S::StateLHA, var::VariableAutomaton) = (S.l_var)[(S.A).map_var_automaton_idx[var]]
+setindex!(S::StateLHA, val::Float64, var::VariableAutomaton) = (S.l_var)[(S.A).map_var_automaton_idx[var]] = val
+setindex!(S::StateLHA, val::Int, var::VariableAutomaton) = (S.l_var)[(S.A).map_var_automaton_idx[var]] = convert(Float64, val)
+function Base.show(io::IO, S::StateLHA)
+ print(io, "State of LHA\n")
+ print(io, "- location: $(S.loc)\n")
+ print(io, "- time: $(S.time)\n")
+ print(io, "- variables:\n")
+ for (var, idx) in (S.A).map_var_automaton_idx
+ print(io, "* $var = $(S.l_var[idx])\n")
+ end
+end
+
+function init_state(A::LHA, x0::SubArray{Int,1}, t0::Float64)
+ S0 = StateLHA(A, "", zeros(length_var(A)), t0)
+ for loc in A.l_loc_init
+ if A.Λ[loc](A,S0)
+ S0.loc = loc
+ break
+ end
+ end
+ return S0
+end
+
+function next_state!(Snplus1::StateLHA, A::LHA,
+ xnplus1::SubArray{Int,1}, tnplus1::Float64, tr_nplus1::Transition,
+ Sn::StateLHA; verbose = false)
+ for i in eachindex(Snplus1.l_var)
+ Snplus1.l_var[i] += (A.l_flow[Sn.loc])[i]*(tnplus1 - Sn.time)
+ end
+ Snplus1.time = tnplus1
+
+ # En fait d'apres observation de Cosmos, après qu'on ait lu la transition on devrait stop.
+ edge_candidates = Edge[]
+ tuple_candidates = Tuple{Location, Location}[]
+ first_round = true
+ detected_event = (tr_nplus1 == nothing) ? true : false
+ turns = 1
+ while first_round || !detected_event || length(edge_candidates) > 0
+ edge_candidates = Edge[]
+ tuple_candidates = Tuple{Location,Location}[]
+ current_loc = Snplus1.loc
+ if verbose
+ @show turns
+ end
+ # Save all edges that satisfies transition predicate (synchronous or autonomous)
+ for loc in A.l_loc
+ tuple_edges = (current_loc, loc)
+ if haskey(A.map_edges, tuple_edges)
+ for edge in A.map_edges[tuple_edges]
+ if edge.check_constraints(A, Snplus1)
+ push!(edge_candidates, edge)
+ push!(tuple_candidates, tuple_edges)
+ end
+ end
+ end
+ end
+ # Search the one we must chose
+ ind_edge = 0
+ for (i, edge) in enumerate(edge_candidates)
+ if edge.transitions[1] == nothing
+ ind_edge = i
+ break
+ end
+ if first_round || !detected_event
+ if (length(edge.transitions) == 1 && tr_nplus1 != nothing && edge.transitions[1] == "ALL") ||
+ (tr_nplus1 in edge.transitions)
+ detected_event = true
+ ind_edge = i
+ end
+ end
+ end
+ # Update the state with the chosen one (if it exists)
+ if ind_edge > 0
+ edge_candidates[ind_edge].update_state!(A, Snplus1, xnplus1)
+ # Should add something like if edges_candidates[ind_edge].transition != nohting break end ??
+ end
+ if verbose
+ @show first_round, detected_event
+ @show tnplus1, tr_nplus1, xnplus1
+ @show ind_edge
+ @show edge_candidates
+ @show tuple_candidates
+ @show Snplus1
+ end
+ if ind_edge == 0 break end
+ # For debug
+ if turns > 10
+ println("Turns, Bad behavior of region2 automaton")
+ @show first_round, detected_event
+ @show length(edge_candidates)
+ @show tnplus1, tr_nplus1, xnplus1
+ @show edge_candidates
+ for edge in edge_candidates
+ @show edge.check_constraints(A, Snplus1)
+ end
+ error("Unpredicted behavior automaton F v2")
+ end
+ turns += 1
+ first_round = false
+ end
+end
+
+# For debug purposes
+function read_trajectory(A::LHA, σ::Trajectory; verbose = false)
+ A_new = LHA(A, (σ.m)._map_obs_var_idx)
+ l_t = times(σ)
+ l_tr = transitions(σ)
+ Sn = init_state(A_new, σ[1], l_t[1])
+ Snplus1 = copy(Sn)
+ if verbose println("Init: ") end
+ if verbose @show Sn end
+ for n in 1:length_states(σ)
+ next_state!(Snplus1, A_new, σ[n], l_t[n], l_tr[n], Sn; verbose = verbose)
+ if Snplus1.loc in A_new.l_loc_final
+ break
+ end
+ Sn = Snplus1
+ end
+ return Sn
+end
+
diff --git a/core/model.jl b/core/model.jl
index e2b23b0..6b46148 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -3,14 +3,16 @@ import StaticArrays: SVector
abstract type Model end
abstract type DiscreteTimeModel <: Model end
+const Transition = Union{String,Nothing}
+# Model types
mutable struct ContinuousTimeModel <: Model
d::Int # state space dim
k::Int # parameter space dim
map_var_idx::Dict{String,Int} # maps str to full state space
_map_obs_var_idx::Dict{String,Int} # maps str to observed state space
map_param_idx::Dict{String,Int} # maps str in parameter space
- l_name_transitions::Vector{String}
+ l_transitions::Vector{Transition}
p::Vector{Float64}
x0::Vector{Int}
t0::Float64
@@ -22,7 +24,7 @@ mutable struct ContinuousTimeModel <: Model
buffer_size::Int
end
-function ContinuousTimeModel(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_name_transitions::Vector{String},
+function ContinuousTimeModel(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_transitions::Vector{String},
p::Vector{Float64}, x0::Vector{Int}, t0::Float64,
f!::Function, is_absorbing::Function;
g::Vector{String} = keys(map_var_idx), time_bound::Float64 = Inf, buffer_size::Int = 10)
@@ -41,9 +43,17 @@ function ContinuousTimeModel(d::Int, k::Int, map_var_idx::Dict, map_param_idx::D
@warn "You have possibly redefined a function Model.is_absorbing used in a previously instantiated model."
end
- return ContinuousTimeModel(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_name_transitions, p, x0, t0, f!, g, _g_idx, is_absorbing, time_bound, buffer_size)
+ return ContinuousTimeModel(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_transitions, p, x0, t0, f!, g, _g_idx, is_absorbing, time_bound, buffer_size)
end
+#=
+mutable struct SynchronizedModel
+ m::ContinuousTimeModel
+ automaton::LHA
+end
+=#
+
+# Simulation
function simulate(m::ContinuousTimeModel)
# trajectory fields
full_values = Matrix{Int}(undef, 1, m.d)
@@ -89,6 +99,62 @@ function simulate(m::ContinuousTimeModel)
return Trajectory(m, values, times, transitions)
end
+#=
+function simulate(product::SynchronizedModel)
+ # trajectory fields
+ m = product.m
+ A = product.automaton
+ full_values = Matrix{Int}(undef, 1, m.d)
+ full_values[1,:] = m.x0
+ times = Float64[m.t0]
+ transitions = Union{String,Nothing}[nothing]
+ reshaped_x0 = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
+ S0 = init_state(A, reshaped_x0, t0)
+ # values at time n
+ n = 0
+ xn = reshaped_x0
+ tn = m.t0
+ Sn = copy(S0)
+ # at time n+1
+ mat_x = zeros(Int, m.buffer_size, m.d)
+ l_t = zeros(Float64, m.buffer_size)
+ l_tr = Vector{String}(undef, m.buffer_size)
+ is_absorbing::Bool = m.is_absorbing(m.p,xn)
+ Snplus1 = copy(Sn)
+ while !is_absorbing && (tn < m.time_bound)
+ i = 0
+ while i < m.buffer_size && !is_absorbing && (tn < m.time_bound)
+ i += 1
+ m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
+ xn = view(mat_x, i, :)
+ tn = l_t[i]
+ tr_n = l_tr[n]
+ A.next_state!(Snplus1, A, xn, tn, tr_n, Sn)
+ Sn = Snplus1
+ is_absorbing = m.is_absorbing(m.p,xn)
+ end
+ full_values = vcat(full_values, view(mat_x, 1:i, :))
+ append!(times, view(l_t, 1:i))
+ append!(transitions, view(l_tr, 1:i))
+ n += i
+ is_absorbing = m.is_absorbing(m.p,xn)
+ end
+ if is_bounded(m)
+ if times[end] > m.time_bound
+ full_values[end,:] = full_values[end-1,:]
+ times[end] = m.time_bound
+ transitions[end] = nothing
+ else
+ full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d))
+ push!(times, m.time_bound)
+ push!(transitions, nothing)
+ end
+ end
+ values = view(full_values, :, m._g_idx)
+ return Trajectory(m, values, times, transitions)
+end
+=#
+
function simulate(m::ContinuousTimeModel, n::Int)
obs = ContinuousObservations(undef, n)
for i = 1:n
diff --git a/core/trajectory.jl b/core/trajectory.jl
index 6f8ce4a..37e9697 100644
--- a/core/trajectory.jl
+++ b/core/trajectory.jl
@@ -6,7 +6,7 @@ struct Trajectory <: AbstractTrajectory
m::ContinuousTimeModel
values::Matrix{Int}
times::Vector{Float64}
- transitions::Vector{Union{String,Nothing}}
+ transitions::Vector{Transition}
end
# Operations
@@ -15,7 +15,7 @@ function -(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end
# Top-level Lp distance function
"""
- `dist_lp(l_σ1, l_σ2; verbose, p, str_stat_list, str_stat_trajectory)`
+`dist_lp(l_σ1, l_σ2; verbose, p, str_stat_list, str_stat_trajectory)`
Function that computes Lp distance between two set of any dimensional trajectories.
...
@@ -35,7 +35,7 @@ function dist_lp(l_σ1::Vector{<:AbstractTrajectory}, l_σ2::Vector{<:AbstractTr
end
"""
- `dist_lp(σ1, σ2; verbose, p, str_stat)`
+`dist_lp(σ1, σ2; verbose, p, str_stat)`
Function that computes Lp distance between two trajectories of any dimension.
It computes Lp distances for each observed variable (contained in `get_obs_var(σ)`). and then apply a statistic on these distances.
@@ -57,7 +57,7 @@ function dist_lp(σ1::AbstractTrajectory, σ2::AbstractTrajectory;
end
"""
- `dist_lp(σ1, σ2, var; verbose, p, str_stat)`
+`dist_lp(σ1, σ2, var; verbose, p, str_stat)`
Function that computes Lp distance between two trajectories of any dimension.
It computes Lp distances for each observed variable (contained in `get_obs_var(σ)`). and then apply a statistic on these distances.
@@ -131,7 +131,7 @@ function dist_lp_mean(σ1::AbstractTrajectory, σ2::AbstractTrajectory;
end
# For a list of trajectories
function dist_lp_mean(l_σ1::Vector{<:AbstractTrajectory}, l_σ2::Vector{<:AbstractTrajectory};
- verbose::Bool = false, p::Int = 1, str_stat_trajectory::String = "mean")
+ verbose::Bool = false, p::Int = 1, str_stat_trajectory::String = "mean")
res = 0.0
for σ1 in l_σ1
for σ2 in l_σ2
@@ -164,12 +164,12 @@ function _riemann_sum(f::Function, t_begin::Real, t_end::Real, step::Float64)
end
function check_consistency(σ::AbstractTrajectory)
- test_sizes = length(σ.times) == length(σ.transitions) == size(σ.values)[1]
- test_obs_var = length_obs_var(σ) == length(σ.m.g) == size(σ.values)[2]
- if !test_sizes error("Issue between sizes of values/times/transitions in this trajectory.") end
- if !test_obs_var error("Issue between sizes of observed variables in model and values") end
+ @assert length(σ.times) == length(σ.transitions) == size(σ.values)[1]
+ @assert length_obs_var(σ) == length(σ.m.g) == size(σ.values)[2]
return true
end
+is_steadystate(σ::AbstractTrajectory) = (σ.m).is_absorbing((σ.m).p, σ[end])
+
# Properties of the trajectory
length_states(σ::AbstractTrajectory) = length(σ.times)
@@ -178,15 +178,29 @@ get_obs_var(σ::AbstractTrajectory) = (σ.m).g
is_bounded(σ::AbstractTrajectory) = σ.transitions[end] == nothing
# Access to trajectory values
-get_var_values(σ::AbstractTrajectory, var::String) =
- view(σ.values, :, (σ.m)._map_obs_var_idx[var])
-get_state(σ::AbstractTrajectory, idx::Int) =
- view(σ.values, idx, :)
-get_value(σ::AbstractTrajectory, var::String, idx::Int) =
- σ.values[idx,(σ.m)._map_obs_var_idx[var]]
-
-
-δ(σ::AbstractTrajectory,t::Float64) = true
+get_var_values(σ::AbstractTrajectory, var::String) = view(σ.values, :, (σ.m)._map_obs_var_idx[var])
+get_state(σ::AbstractTrajectory, idx::Int) = view(σ.values, idx, :)
+get_value(σ::AbstractTrajectory, var::String, idx::Int) = σ.values[idx,(σ.m)._map_obs_var_idx[var]]
+# Operation @
+function get_state_from_time(σ::AbstractTrajectory, t::Float64)
+ @assert t >= 0.0
+ l_t = times(σ)
+ if t == l_t[end] return σ[end] end
+ if t > l_t[end]
+ if !is_bounded(σ)
+ return σ[end]
+ else
+ error("This trajectory is bounded and you're accessing out of the bounds")
+ end
+ end
+ for i in eachindex(l_t)
+ if l_t[i] <= t < l_t[i+1]
+ return σ[i]
+ end
+ end
+ error("Unexpected behavior")
+end
+δ(σ::AbstractTrajectory,idx::Int) = times(σ)[i+1] - times(σ)[i]
states(σ::AbstractTrajectory) = σ.values
times(σ::AbstractTrajectory) = σ.times
@@ -198,4 +212,5 @@ getindex(σ::AbstractTrajectory, idx::Int) = get_state(σ, idx)
getindex(σ::AbstractTrajectory, var::String, idx::Int) = get_value(σ, var, idx)
# Get the path of a variable ["I"]
getindex(σ::AbstractTrajectory, var::String) = get_var_values(σ, var)
+lastindex(σ::AbstractTrajectory) = length_states(σ)
diff --git a/tests/automata/automaton_F.jl b/tests/automata/automaton_F.jl
new file mode 100644
index 0000000..8b70061
--- /dev/null
+++ b/tests/automata/automaton_F.jl
@@ -0,0 +1,12 @@
+
+using MarkovProcesses
+
+load_model("SIR")
+load_automaton("automaton_F")
+SIR.time_bound = 120.0
+x1, x2, t1, t2 = 0.0, Inf, 100.0, 120.0
+A_F = create_automaton_F(SIR, x1, x2, t1, t2, "I") # <: LHA
+
+#sync_SIR = SIR * A_F
+#σ, state_lha = simulate(sync_SIR)
+
diff --git a/tests/automata/read_trajectory_last_state_F.jl b/tests/automata/read_trajectory_last_state_F.jl
new file mode 100644
index 0000000..636080d
--- /dev/null
+++ b/tests/automata/read_trajectory_last_state_F.jl
@@ -0,0 +1,32 @@
+
+using MarkovProcesses
+
+load_model("SIR")
+load_automaton("automaton_F")
+SIR.time_bound = 120.0
+x1, x2, t1, t2 = 0.0, Inf, 100.0, 120.0
+
+A_F = create_automaton_F(SIR, x1, x2, t1, t2, "I") # <: LHA
+
+function test_last_state(A::LHA, m::ContinuousTimeModel)
+ σ = simulate(m)
+ Send = read_trajectory(A, σ)
+ test = (get_state_from_time(σ, Send.time)[1] == Send["n"]) & (Send["d"] == 0)
+ #=
+ if !test
+ @show Send
+ @show get_state_from_time(σ, Send.time)
+ error("tkt")
+ end=#
+ return test
+end
+
+test_all = true
+nbr_sim = 10000
+for i = 1:nbr_sim
+ test = test_last_state(A_F, SIR)
+ global test_all = test_all & test
+end
+
+return test_all
+
diff --git a/tests/run_all.jl b/tests/run_all.jl
index e1433d2..5ed6f29 100644
--- a/tests/run_all.jl
+++ b/tests/run_all.jl
@@ -2,4 +2,5 @@
include("run_unit.jl")
include("run_simulation.jl")
include("run_dist_lp.jl")
+include("run_automata.jl")
diff --git a/tests/run_automata.jl b/tests/run_automata.jl
new file mode 100644
index 0000000..abec4d0
--- /dev/null
+++ b/tests/run_automata.jl
@@ -0,0 +1,7 @@
+
+using Test
+
+@testset "Automata tests" begin
+ @test include("automata/read_trajectory_last_state_F.jl")
+end
+
diff --git a/tests/run_bench.jl b/tests/run_bench.jl
new file mode 100644
index 0000000..5b6a694
--- /dev/null
+++ b/tests/run_bench.jl
@@ -0,0 +1,16 @@
+
+using Test
+import MarkovProcesses: get_module_path
+
+path_bench = get_module_path() * "/bench/"
+function include_arg!(path_file::String, arg::String)
+ global ARGS = [arg]
+ include(path_bench * "array_memory_order/sim.jl")
+ return true
+end
+
+@testset "Benchmark" begin
+ @test include_arg!(path_bench * "array_memory_order/sim.jl", "SIR")
+ @test include_arg!(path_bench * "pygmalion/sim_sir.jl", "SIR")
+end
+
--
GitLab