diff --git a/automata/automaton_F.jl b/automata/automaton_F.jl
new file mode 100644
index 0000000000000000000000000000000000000000..38c2c3e3ec5d91958ebb53f7a1292fa58a98ec90
--- /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 e9d26143b64737cb6bfd92bb7230dffc295c5ba6..dfb23e636f9006dae1ff89f483f3f3068cf29155 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 0000000000000000000000000000000000000000..c3a1359f21325003170185931e2a73f370e35d9a
--- /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 e2b23b07b0b24c4bac76e3449a3cc51f3a47b75d..6b4614814c5ba92a00634483ede1e95337fc5116 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 6f8ce4ad82bd4a37f03a74d3ef35e54a6b54ef87..37e969764f16ad6281e3ded37f380b6a5076c510 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 0000000000000000000000000000000000000000..8b70061e297bb87ed66ce9027926ee051b1ed267
--- /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 0000000000000000000000000000000000000000..636080d5abe3e04ccbc4ad9b18f1ab115c47cede
--- /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 e1433d2a44cd2e7de93f15912c657758868ad269..5ed6f292c64e5b11d3a0008f10150b3ec9ec2cf0 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 0000000000000000000000000000000000000000..abec4d04c8ce02da9f18e5067be5dc6eceaf93d9
--- /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 0000000000000000000000000000000000000000..5b6a694fa7e58d57e50e3f5d8f0a80e9b78b4d1d
--- /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
+