Commit 21c72d7e authored by Bentriou Mahmoud's avatar Bentriou Mahmoud
Browse files

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.
parent cefe75b7
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
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
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
......@@ -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
......
......@@ -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(σ)
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)
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
......@@ -2,4 +2,5 @@
include("run_unit.jl")
include("run_simulation.jl")
include("run_dist_lp.jl")
include("run_automata.jl")
using Test
@testset "Automata tests" begin
@test include("automata/read_trajectory_last_state_F.jl")
end
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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment