From f4c41acb822d97b116aeffc9efe5923e9d5e3e72 Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Tue, 1 Dec 2020 18:04:17 +0100
Subject: [PATCH] Test of distance G automaton passes!
The two Cosmos tests checks the statistical correctness of estimated
distance value for R1,R2,R3,R5 experiments. Also, they test that all
trajectories are accepted.
A bunch of errors were fixed in automata files (Cosmos .lha and julia
./automata files) and simulation function that introduced statistical
bias, i.e. approximated values that were closed but not in the
confidence interval. The general structure of simulation and was
reworked in order to suits well te behavior of Cosmos.
---
automata/automaton_F.jl | 57 ++++++++----
automata/automaton_G.jl | 94 +++++++++-----------
core/lha.jl | 122 +++++++++++++++++++-------
core/model.jl | 53 +++++++----
core/plots.jl | 3 +-
tests/cosmos/distance_F/ER_1D.jl | 6 +-
tests/cosmos/distance_F/dist_F_ER.lha | 23 +++--
tests/cosmos/distance_G/ER_R5.jl | 12 ++-
tests/cosmos/distance_G/dist_G_ER.lha | 2 +-
tests/run_cosmos.jl | 1 +
10 files changed, 239 insertions(+), 134 deletions(-)
diff --git a/automata/automaton_F.jl b/automata/automaton_F.jl
index b464a7a..1918388 100644
--- a/automata/automaton_F.jl
+++ b/automata/automaton_F.jl
@@ -35,43 +35,62 @@ function create_automaton_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
tuple = ("l0", "l1")
cc_aut_F_l0l1_1(A::LHA, S::StateLHA) = true
us_aut_F_l0l1_1!(A::LHA, S::StateLHA, x::Vector{Int}) =
- (S.loc = "l1"; S["d"] = Inf; S["n"] = get_value(A, x, str_obs))
+ (S.loc = "l1";
+ S["d"] = Inf;
+ S["n"] = get_value(A, x, str_obs);
+ S["isabs"] = m.isabsorbing(m.p,x))
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::Vector{Int}) = (S["d"] = 0; S.loc = "l2")
+ (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::Vector{Int}) =
+ (S.loc = "l2";
+ S["d"] = 0)
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"] ||Â istrue(S["isabs"])))
- us_aut_F_l1l2_2!(A::LHA, S::StateLHA, x::Vector{Int}) = (S.loc = "l2")
+ cc_aut_F_l1l2_2(A::LHA, S::StateLHA) =
+ S["d"] > 0 &&
+ (S.time > A.l_ctes["t2"] ||Â istrue(S["isabs"]))
+ us_aut_F_l1l2_2!(A::LHA, S::StateLHA, x::Vector{Int}) =
+ (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::Vector{Int}) = (S.loc = "l2")
+ 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::Vector{Int}) =
+ (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::Vector{Int}) = (S["d"] = 0; S.loc = "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::Vector{Int}) =
+ (S.loc = "l3";
+ S["d"] = 0;)
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"])
+ (S["n"] < A.l_ctes["x1"] || S["n"] > A.l_ctes["x2"]) &&
+ (S.time <= A.l_ctes["t1"])
us_aut_F_l1l3_2!(A::LHA, S::StateLHA, x::Vector{Int}) =
- (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")
+ (S.loc = "l3";
+ 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)))
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"])
+ (S["n"] < A.l_ctes["x1"] || 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::Vector{Int}) =
- (S["d"] = min(S["d"], min(abs(S["n"] - A.l_ctes["x1"]), abs(S["n"] - A.l_ctes["x2"]))); S.loc = "l3")
+ (S.loc = "l3";
+ S["d"] = min(S["d"], min(abs(S["n"] - A.l_ctes["x1"]), abs(S["n"] - A.l_ctes["x2"]))))
edge3 = Edge([nothing], cc_aut_F_l1l3_3, us_aut_F_l1l3_3!)
map_edges[tuple] = [edge1, edge2, edge3]
@@ -79,11 +98,15 @@ function create_automaton_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
tuple = ("l3", "l1")
cc_aut_F_l3l1_1(A::LHA, S::StateLHA) = true
us_aut_F_l3l1_1!(A::LHA, S::StateLHA, x::Vector{Int}) =
- (S["n"] = get_value(A, x, str_obs); S.loc = "l1"; S["isabs"] = m.isabsorbing(m.p,x))
+ (S.loc = "l1";
+ S["n"] = get_value(A, x, str_obs);
+ S["isabs"] = m.isabsorbing(m.p,x))
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::Vector{Int}) = (S["n"] = get_value(A, x, str_obs); S.loc = "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::Vector{Int}) =
+ (S.loc = "l2")
edge2 = Edge([nothing], cc_aut_F_l3l2_1, us_aut_F_l3l2_1!)
map_edges[tuple] = [edge1, edge2]
diff --git a/automata/automaton_G.jl b/automata/automaton_G.jl
index 9c982c9..c221098 100644
--- a/automata/automaton_G.jl
+++ b/automata/automaton_G.jl
@@ -15,7 +15,7 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
l_loc_final = ["l2"]
## Map of automaton variables
- map_var_automaton_idx = Dict{VariableAutomaton,Int}("t'" => 1, "in" => 2,
+ map_var_automaton_idx = Dict{VariableAutomaton,Int}("tprime" => 1, "in" => 2,
"n" => 3, "d" => 4, "isabs" => 5)
## Flow of variables
@@ -34,15 +34,19 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
tuple = ("l0", "l1")
cc_aut_G_l0l1_1(A::LHA, S::StateLHA) = true
us_aut_G_l0l1_1!(A::LHA, S::StateLHA, x::Vector{Int}) =
- (S.loc = "l1"; S["d"] = 0; S["n"] = get_value(A, x, str_obs); S["in"] = true; S["isabs"] = m.isabsorbing(m.p,x))
+ (S.loc = "l1";
+ S["d"] = 0;
+ S["n"] = get_value(A, x, str_obs);
+ S["in"] = true;
+ S["isabs"] = m.isabsorbing(m.p,x))
edge1 = Edge([nothing], cc_aut_G_l0l1_1, us_aut_G_l0l1_1!)
map_edges[tuple] = [edge1]
# l1 loc
tuple = ("l1", "l3")
cc_aut_G_l1l3_1(A::LHA, S::StateLHA) =
- (S.time < A.l_ctes["t1"] &&
- (S["n"] < A.l_ctes["x1"] || S["n"] > A.l_ctes["x2"]))
+ S.time <= A.l_ctes["t1"] &&
+ S["n"] < A.l_ctes["x1"] || S["n"] > A.l_ctes["x2"]
us_aut_G_l1l3_1!(A::LHA, S::StateLHA, x::Vector{Int}) =
(S.loc = "l3";
S["d"] = min(abs(A.l_ctes["x1"] - S["n"]), abs(A.l_ctes["x2"] - S["n"]));
@@ -56,7 +60,7 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
us_aut_G_l1l3_3!(A::LHA, S::StateLHA, x::Vector{Int}) =
(S.loc = "l3";
S["d"] = S["d"] * (S.time - A.l_ctes["t1"]);
- S["t'"] = 0.0)
+ S["tprime"] = 0.0)
edge3 = Edge([nothing], cc_aut_G_l1l3_3, us_aut_G_l1l3_3!)
cc_aut_G_l1l3_2(A::LHA, S::StateLHA) =
@@ -74,7 +78,7 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
(A.l_ctes["x1"] <= S["n"] <= A.l_ctes["x2"])
us_aut_G_l1l3_4!(A::LHA, S::StateLHA, x::Vector{Int}) =
(S.loc = "l3";
- S["t'"] = 0.0)
+ S["tprime"] = 0.0)
edge4 = Edge([nothing], cc_aut_G_l1l3_4, us_aut_G_l1l3_4!)
map_edges[tuple] = [edge1, edge2, edge3, edge4]
@@ -111,7 +115,23 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
(S.loc = "l2";
S["d"] = S["d"] * (A.l_ctes["t2"] - A.l_ctes["t1"]))
edge2 = Edge([nothing], cc_aut_G_l1l2_2, us_aut_G_l1l2_2!)
- map_edges[tuple] = [edge1, edge2]
+ cc_aut_G_l1l2_3(A::LHA, S::StateLHA) =
+ istrue(S["isabs"]) &&
+ S.time <= A.l_ctes["t1"]
+ us_aut_G_l1l2_3!(A::LHA, S::StateLHA, x::Vector{Int}) =
+ (S.loc = "l2";
+ S["d"] = S["d"] * (A.l_ctes["t2"] - A.l_ctes["t1"]))
+ edge3 = Edge([nothing], cc_aut_G_l1l2_3, us_aut_G_l1l2_3!)
+ cc_aut_G_l1l2_4(A::LHA, S::StateLHA) =
+ istrue(S["isabs"]) &&
+ (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"])
+ us_aut_G_l1l2_4!(A::LHA, S::StateLHA, x::Vector{Int}) =
+ (S.loc = "l2";
+ S["d"] += (A.l_ctes["t2"] - S.time) *
+ min(abs(A.l_ctes["x1"] - S["n"]), abs(A.l_ctes["x1"] - S["n"])))
+ edge4 = Edge([nothing], cc_aut_G_l1l2_4, us_aut_G_l1l2_4!)
+
+ map_edges[tuple] = [edge1, edge2, edge3, edge4]
# l3 loc
tuple = ("l3", "l1")
@@ -124,44 +144,29 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
map_edges[tuple] = [edge1]
tuple = ("l3", "l2")
- cc_aut_G_l3l2_1(A::LHA, S::StateLHA) =
- istrue(S["in"]) &&
- (S.time >= A.l_ctes["t2"] || istrue(S["isabs"]))
- us_aut_G_l3l2_1!(A::LHA, S::StateLHA, x::Vector{Int}) =
- (S.loc = "l2")
- edge1 = Edge([nothing], cc_aut_G_l3l2_1, us_aut_G_l3l2_1!)
cc_aut_G_l3l2_2(A::LHA, S::StateLHA) =
- !istrue(S["in"]) &&
- (S.time >= A.l_ctes["t2"] || istrue(S["isabs"]))
+ istrue(S["in"]) &&
+ S.time >= A.l_ctes["t2"]
us_aut_G_l3l2_2!(A::LHA, S::StateLHA, x::Vector{Int}) =
(S.loc = "l2";
S["d"] = S["d"] * (A.l_ctes["t2"] - A.l_ctes["t1"]))
edge2 = Edge([nothing], cc_aut_G_l3l2_2, us_aut_G_l3l2_2!)
- cc_aut_G_l3l2_3(A::LHA, S::StateLHA) =
- istrue(S["isabs"]) &&
- S.time <= A.l_ctes["t1"]
- us_aut_G_l3l2_3!(A::LHA, S::StateLHA, x::Vector{Int}) =
- (S.loc = "l2";
- S["d"] = S["d"] * (A.l_ctes["t2"] - A.l_ctes["t1"]))
- edge3 = Edge([nothing], cc_aut_G_l3l2_3, us_aut_G_l3l2_3!)
- cc_aut_G_l3l2_4(A::LHA, S::StateLHA) =
- istrue(S["isabs"]) &&
- S.time >= A.l_ctes["t1"]
- us_aut_G_l3l2_4!(A::LHA, S::StateLHA, x::Vector{Int}) =
- (S.loc = "l2";
- S["d"] += (A.l_ctes["t2"] - A.l_ctes["t1"]) *
- min(abs(S["n"] - A.l_ctes["x1"]), abs(S["n"] - A.l_ctes["x2"])))
- edge4 = Edge([nothing], cc_aut_G_l3l2_4, us_aut_G_l3l2_4!)
+ cc_aut_G_l3l2_1(A::LHA, S::StateLHA) =
+ !istrue(S["in"]) &&
+ S.time >= A.l_ctes["t2"]
+ us_aut_G_l3l2_1!(A::LHA, S::StateLHA, x::Vector{Int}) =
+ (S.loc = "l2")
+ edge1 = Edge([nothing], cc_aut_G_l3l2_1, us_aut_G_l3l2_1!)
- map_edges[tuple] = [edge1, edge2, edge3, edge4]
+ map_edges[tuple] = [edge1, edge2]
# l4 loc
tuple = ("l4", "l1")
cc_aut_G_l4l1_1(A::LHA, S::StateLHA) = true
us_aut_G_l4l1_1!(A::LHA, S::StateLHA, x::Vector{Int}) =
(S.loc = "l1";
- S["d"] += S["t'"] * min(abs(A.l_ctes["x1"] - S["n"]), abs(A.l_ctes["x2"] - S["n"]));
- S["t'"] = 0.0;
+ S["d"] += S["tprime"] * min(abs(A.l_ctes["x1"] - S["n"]), abs(A.l_ctes["x2"] - S["n"]));
+ S["tprime"] = 0.0;
S["n"] = get_value(A, x, str_obs);
S["in"] = true;
S["isabs"] = m.isabsorbing(m.p,x))
@@ -170,29 +175,14 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
tuple = ("l4", "l2")
cc_aut_G_l4l2_1(A::LHA, S::StateLHA) =
- (S.time >= A.l_ctes["t2"] || istrue(S["isabs"]))
+ (S.time >= A.l_ctes["t2"])
us_aut_G_l4l2_1!(A::LHA, S::StateLHA, x::Vector{Int}) =
(S.loc = "l2";
- S["d"] += S["t'"] * min(abs(A.l_ctes["x1"] - S["n"]), abs(A.l_ctes["x2"] - S["n"]));
- S["t'"] = 0.0)
+ S["d"] += S["tprime"] * min(abs(A.l_ctes["x1"] - S["n"]), abs(A.l_ctes["x2"] - S["n"]));
+ S["tprime"] = 0.0)
edge1 = Edge([nothing], cc_aut_G_l4l2_1, us_aut_G_l4l2_1!)
- cc_aut_G_l4l2_2(A::LHA, S::StateLHA) =
- istrue(S["isabs"]) &&
- S.time <= A.l_ctes["t1"]
- us_aut_G_l4l2_2!(A::LHA, S::StateLHA, x::Vector{Int}) =
- (S.loc = "l2";
- S["d"] = S["d"] * (A.l_ctes["t2"] - A.l_ctes["t1"]))
- edge2 = Edge([nothing], cc_aut_G_l4l2_2, us_aut_G_l4l2_2!)
- cc_aut_G_l4l2_3(A::LHA, S::StateLHA) =
- istrue(S["isabs"]) &&
- S.time >= A.l_ctes["t1"]
- us_aut_G_l4l2_3!(A::LHA, S::StateLHA, x::Vector{Int}) =
- (S.loc = "l2";
- S["d"] += (A.l_ctes["t2"] - A.l_ctes["t1"]) *
- min(abs(S["n"] - A.l_ctes["x1"]), abs(S["n"] - A.l_ctes["x2"])))
- edge3 = Edge([nothing], cc_aut_G_l4l2_3, us_aut_G_l4l2_3!)
- map_edges[tuple] = [edge1, edge2, edge3]
+ map_edges[tuple] = [edge1]
## Constants
l_ctes = Dict{String,Float64}("x1" => x1, "x2" => x2, "t1" => t1, "t2" => t2)
diff --git a/core/lha.jl b/core/lha.jl
index fcc3fce..af93550 100644
--- a/core/lha.jl
+++ b/core/lha.jl
@@ -43,34 +43,39 @@ function init_state(A::LHA, x0::Vector{Int}, t0::Float64)
end
function _find_edge_candidates!(edge_candidates::Vector{Edge}, current_loc::Location,
- A::LHA, Snplus1::StateLHA)
-
+ A::LHA, Snplus1::StateLHA, only_asynchronous::Bool)
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)
+ if edge.transitions[1] == nothing
+ pushfirst!(edge_candidates, edge)
+ end
+ if !only_asynchronous && edge.transitions[1] != nothing
+ push!(edge_candidates, edge)
+ end
end
end
end
end
end
-function _get_edge_index(edge_candidates::Vector{Edge}, first_round::Bool,
+function _get_edge_index(edge_candidates::Vector{Edge},
detected_event::Bool, tr_nplus1::Transition)
ind_edge = 0
bool_event = detected_event
for i in eachindex(edge_candidates)
edge = edge_candidates[i]
+ # Asynchronous detection
if edge.transitions[1] == nothing
return (i, bool_event)
end
- if first_round || !detected_event
- if (length(edge.transitions) == 1 && tr_nplus1 != nothing && edge.transitions[1] == "ALL") ||
- (tr_nplus1 in edge.transitions)
- ind_edge = i
- bool_event = true
+ # Synchronous detection
+ if !detected_event && tr_nplus1 != nothing
+ if (length(edge.transitions) == 1 && edge.transitions[1] == "ALL") ||
+ (tr_nplus1 in edge.transitions)
+ return (i, true)
end
end
end
@@ -80,41 +85,45 @@ end
function next_state!(Snplus1::StateLHA, A::LHA,
xnplus1::Vector{Int}, tnplus1::Float64, tr_nplus1::Transition,
Sn::StateLHA; verbose::Bool = 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[]
first_round::Bool = true
detected_event::Bool = (tr_nplus1 == nothing) ? true : false
- turns = 1
- while first_round || !detected_event || length(edge_candidates) > 0
+ turns = 0
+
+ if verbose
+ println("#####")
+ @show xnplus1, tnplus1, tr_nplus1
+ @show Sn
+ @show Snplus1
+ end
+
+ # In terms of values not reference, Snplus1 == Sn
+ # First, we check the asynchronous transitions
+ while first_round || length(edge_candidates) > 0
+ turns += 1
edge_candidates = empty!(edge_candidates)
current_loc = Snplus1.loc
- if verbose @show turns end
- # Save all edges that satisfies transition predicate (synchronous or autonomous)
- _find_edge_candidates!(edge_candidates, current_loc, A, Snplus1)
- # Search the one we must chose
- ind_edge, detected_event = _get_edge_index(edge_candidates, first_round, detected_event, tr_nplus1)
+ # Save all edges that satisfies transition predicate (asynchronous ones)
+ _find_edge_candidates!(edge_candidates, current_loc, A, Snplus1, 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, detected_event, nothing)
# 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 (ind_edge == 0)
- #if (ind_edge == 0 ||Â detected_event)
- break
end
+ first_round = false
if verbose
- @show first_round detected_event
- @show tnplus1 tr_nplus1 xnplus1
- @show ind_edge
+ @show turns
@show edge_candidates
- @show tuple_candidates
+ @show ind_edge, detected_event
+ println("After update")
@show Snplus1
end
+ if (ind_edge == 0)
+ break
+ end
# For debug
if turns > 100
println("Number of turns in next_state! is suspicious")
@@ -127,8 +136,61 @@ function next_state!(Snplus1::StateLHA, A::LHA,
end
error("Unpredicted behavior automaton")
end
+ end
+ if verbose
+ @show Snplus1
+ println("Time flies with the flow...")
+ end
+ # Now time flies according to the flow
+ for i in eachindex(Snplus1.l_var)
+ coeff_deriv = (A.l_flow[Snplus1.loc])[i]
+ if coeff_deriv > 0
+ Snplus1.l_var[i] += coeff_deriv*(tnplus1 - Snplus1.time)
+ end
+ end
+ Snplus1.time = tnplus1
+ if verbose
+ @show Snplus1
+ end
+ # Now firing an edge according to the event
+ first_round = true
+ while first_round || length(edge_candidates) > 0
turns += 1
+ edge_candidates = empty!(edge_candidates)
+ current_loc = Snplus1.loc
+ # Save all edges that satisfies transition predicate (synchronous ones)
+ _find_edge_candidates!(edge_candidates, current_loc, A, Snplus1, false)
+ # Search the one we must chose
+ ind_edge, detected_event = _get_edge_index(edge_candidates, detected_event, tr_nplus1)
+ # Update the state with the chosen one (if it exists)
+ if verbose
+ @show turns
+ @show edge_candidates
+ @show ind_edge, detected_event
+ end
+ if ind_edge > 0
+ edge_candidates[ind_edge].update_state!(A, Snplus1, xnplus1)
+ end
first_round = false
+ if verbose
+ println("After update")
+ @show Snplus1
+ end
+ if (ind_edge == 0 ||Â detected_event)
+ break
+ end
+ # For debug
+ if turns > 100
+ println("Number of turns in next_state! is suspicious")
+ @show 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")
+ end
end
end
diff --git a/core/model.jl b/core/model.jl
index da8faf3..c8404b7 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -125,7 +125,8 @@ function simulate(m::ContinuousTimeModel; p::Union{Nothing,AbstractVector{Float6
return Trajectory(m, values, times, transitions)
end
-function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Float64}} = nothing)
+function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Float64}} = nothing,
+ verbose::Bool = false)
m = product.m
A = product.automaton
p_sim = m.p
@@ -149,6 +150,11 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
Sn = copy(S0)
isabsorbing::Bool = m.isabsorbing(p_sim,xn)
isacceptedLHA::Bool = isaccepted(Sn)
+ # Alloc of vectors where we stock n+1 values
+ vec_x = zeros(Int, m.d)
+ l_t = Float64[0.0]
+ l_tr = Transition[nothing]
+ Snplus1 = copy(Sn)
# If x0 is absorbing
if isabsorbing || isacceptedLHA
_resize_trajectory!(full_values, times, transitions, 1)
@@ -156,13 +162,12 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
if isbounded(m)
_finish_bounded_trajectory!(values, times, transitions, m.time_bound)
end
+ if isabsorbing && !isacceptedLHA
+ next_state!(Snplus1, A, xn, m.time_bound, nothing, Sn; verbose = verbose)
+ copyto!(Sn, Snplus1)
+ end
return SynchronizedTrajectory(Sn, product, values, times, transitions)
end
- # Alloc of vectors where we stock n+1 values
- vec_x = zeros(Int, m.d)
- l_t = Float64[0.0]
- l_tr = Transition[nothing]
- Snplus1 = copy(Sn)
# First we fill the allocated array
for i = 2:m.estim_min_states
m.f!(vec_x, l_t, l_tr, xn, tn, p_sim)
@@ -173,7 +178,7 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
n += 1
copyto!(xn, vec_x)
tr_n = l_tr[1]
- next_state!(Snplus1, A, xn, tn, tr_n, Sn)
+ next_state!(Snplus1, A, xn, tn, tr_n, Sn; verbose = verbose)
_update_values!(full_values, times, transitions, xn, tn, tr_n, i)
copyto!(Sn, Snplus1)
isabsorbing = m.isabsorbing(p_sim,xn)
@@ -189,6 +194,10 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
if isbounded(m)
_finish_bounded_trajectory!(values, times, transitions, m.time_bound)
end
+ if isabsorbing && !isacceptedLHA
+ next_state!(Snplus1, A, xn, m.time_bound, nothing, Sn; verbose = verbose)
+ copyto!(Sn, Snplus1)
+ end
return SynchronizedTrajectory(Sn, product, values, times, transitions)
end
# Otherwise, buffering system
@@ -207,7 +216,7 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
end
copyto!(xn, vec_x)
tr_n = l_tr[1]
- next_state!(Snplus1, A, xn, tn, tr_n, Sn)
+ next_state!(Snplus1, A, xn, tn, tr_n, Sn; verbose = verbose)
_update_values!(full_values, times, transitions,
xn, tn, tr_n, m.estim_min_states+size_tmp+i)
copyto!(Sn, Snplus1)
@@ -230,6 +239,10 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
# the trajectory is bounded
_finish_bounded_trajectory!(values, times, transitions, m.time_bound)
end
+ if isabsorbing && !isacceptedLHA
+ next_state!(Snplus1, A, xn, m.time_bound, nothing, Sn; verbose = verbose)
+ copyto!(Sn, Snplus1)
+ end
return SynchronizedTrajectory(Sn, product, values, times, transitions)
end
@@ -249,20 +262,20 @@ function volatile_simulate(product::SynchronizedModel;
Sn = copy(S0)
isabsorbing::Bool = m.isabsorbing(p_sim,xn)
isacceptedLHA::Bool = isaccepted(Sn)
- # If x0 is absorbing
- if isabsorbing || isacceptedLHA
- return Sn
- end
# Alloc of vectors where we stock n+1 values
vec_x = zeros(Int, m.d)
l_t = Float64[0.0]
l_tr = Transition[nothing]
Snplus1 = copy(Sn)
- while !isabsorbing &&Â tn <= m.time_bound && !isacceptedLHA
- if verbose
- @show vec_x
- @show Sn
+ # If x0 is absorbing
+ if isabsorbing || isacceptedLHA
+ if !isacceptedLHA
+ next_state!(Snplus1, A, xn, m.time_bound, nothing, Sn; verbose = verbose)
+ copyto!(Sn, Snplus1)
end
+ return Sn
+ end
+ while !isabsorbing &&Â tn <= m.time_bound && !isacceptedLHA
m.f!(vec_x, l_t, l_tr, xn, tn, p_sim)
tn = l_t[1]
if tn > m.time_bound
@@ -271,12 +284,16 @@ function volatile_simulate(product::SynchronizedModel;
end
copyto!(xn, vec_x)
tr_n = l_tr[1]
- next_state!(Snplus1, A, xn, tn, tr_n, Sn)
+ next_state!(Snplus1, A, xn, tn, tr_n, Sn; verbose = verbose)
copyto!(Sn, Snplus1)
isabsorbing = m.isabsorbing(p_sim,xn)
isacceptedLHA = isaccepted(Snplus1)
n += 1
end
+ if isabsorbing && !isacceptedLHA
+ next_state!(Snplus1, A, xn, m.time_bound, nothing, Sn; verbose = verbose)
+ copyto!(Sn, Snplus1)
+ end
return Sn
end
@@ -313,7 +330,7 @@ function mean_value_lha(sm::SynchronizedModel, str_var::String, nbr_sim::Int)
return sum_val / nbr_sim
end
-function distribute_prob_accept_lha(sm::SynchronizedModel, str_var::String, nbr_sim::Int)
+function distribute_prob_accept_lha(sm::SynchronizedModel, nbr_sim::Int)
sum_val = @distributed (+) for i = 1:nbr_sim Int(isaccepted(volatile_simulate(sm))) end
return sum_val / nbr_sim
end
diff --git a/core/plots.jl b/core/plots.jl
index 9e2ac46..898267c 100644
--- a/core/plots.jl
+++ b/core/plots.jl
@@ -52,8 +52,7 @@ end
function plot!(A::LHA)
x1, x2, t1, t2 = A.l_ctes["x1"], A.l_ctes["x2"], A.l_ctes["t1"], A.l_ctes["t2"]
- rectangle(w, h, x, y) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])
- plot!(rectangle(t2, x2, t1, x1), opacity = 0.5)
+ plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5)
end
export plot, plot!
diff --git a/tests/cosmos/distance_F/ER_1D.jl b/tests/cosmos/distance_F/ER_1D.jl
index 80c8796..3b9176f 100644
--- a/tests/cosmos/distance_F/ER_1D.jl
+++ b/tests/cosmos/distance_F/ER_1D.jl
@@ -42,6 +42,7 @@ for exp in l_exp
--const k_3=$(k3),x1=$x1,x2=$x2,t1=$t1,t2=$t2
--level $(level) --width $(width)
--verbose 0`
+ #run(command)
run(pipeline(command, stderr=devnull))
dict_values = cosmos_get_values("Result_dist_F_$(str_model).res")
l_dist_cosmos[i] = dict_values["Estimated value"]
@@ -52,8 +53,11 @@ for exp in l_exp
set_param!(ER, "k3", convert(Float64, k3))
sync_ER = ER*A_F
l_dist_pkg[i] = distribute_mean_value_lha(sync_ER, "d", nb_sim)
+ nb_accepts_pkg = distribute_prob_accept_lha(sync_ER, nb_sim)
+ #@info "About accepts" nb_sim nb_accepted nb_accepts_pkg
test = isapprox(l_dist_cosmos[i], l_dist_pkg[i]; atol = width*1.01)
- global test_all = test_all && test
+ test2 = nb_accepts_pkg == (nb_sim / nb_accepted)
+ global test_all = test_all && test && test2
if !test
@show l_dist_pkg[i], l_dist_cosmos[i]
@show exp
diff --git a/tests/cosmos/distance_F/dist_F_ER.lha b/tests/cosmos/distance_F/dist_F_ER.lha
index baa9f54..47cd3fb 100644
--- a/tests/cosmos/distance_F/dist_F_ER.lha
+++ b/tests/cosmos/distance_F/dist_F_ER.lha
@@ -1,4 +1,4 @@
-NbVariables = 3;
+NbVariables = 4;
NbLocations = 4;
const int x1 = 50;
@@ -7,7 +7,7 @@ const int x2 = 75;
const double t1 = 0.025;
const double t2 = 0.05;
-VariablesList = {t,d,n};
+VariablesList = {t,d,n,test_abs};
LocationsList = {l0, l1, l2, l3};
@@ -24,17 +24,22 @@ Locations={
};
Edges={
-((l0,l1), #, t>=0, {n=P});
+((l0,l1), #, t>=0, {n=P,test_abs=k_1*(E*S)+k_2*ES});
+
+((l1,l2),#, n>=x1 & n<=x2 & t>=t1 & t<=t2 ,{d=0});
+((l1,l2),#, t>=t2 , #);
+((l1,l2),#, test_abs=0 , #);
+((l1,l2),#, d=0 & t>=t1, #);
+
+((l1,l3), #, n>=x1 & n<=x2, {d=0});
((l1,l3), #, t<=t1 & n<=x1-1, {d=min(((t-t1)^2+(n-x2)^2)^0.5,((t-t1)^2+(n-x1)^2)^0.5)});
((l1,l3), #, t<=t1 & n>=x2+1, {d=min(((t-t1)^2+(n-x2)^2)^0.5,((t-t1)^2+(n-x1)^2)^0.5)});
-((l1,l3), #, n>=x1 & n<=x2, {d=0});
-((l1,l3), #, n<=x1-1 & t >= t1 & t<=t2, {d=min(d ,min( ((n-x1)^2)^0.5, ((n-x2)^2)^0.5 ) )});
+((l1,l3), #, n<=x1-1 & t >= t1 & t<=t2, {d=min(d,min(((n-x1)^2)^0.5,((n-x2)^2)^0.5))});
((l1,l3), #, n>=x2+1 & t >= t1 & t<=t2, {d=min(d,min( ((n-x1)^2)^0.5, ((n-x2)^2)^0.5 ) ) });
-((l3,l1), ALL, t>=0, {n=P});
+
+((l3,l1), ALL, t>=0, {n=P,test_abs=k_1*(E*S)+k_2*ES});
((l3,l2), #, t>=t2 ,#);
-((l1,l2),#, n>=x1 & n<=x2 & t>=t1 & t<=t2 ,{d=0});
-((l1,l2),#, t>=t2 , #);
-((l1,l2),#, d=0 & t>=t1, #);
};
+
diff --git a/tests/cosmos/distance_G/ER_R5.jl b/tests/cosmos/distance_G/ER_R5.jl
index b211a29..c62185c 100644
--- a/tests/cosmos/distance_G/ER_R5.jl
+++ b/tests/cosmos/distance_G/ER_R5.jl
@@ -1,6 +1,7 @@
@everywhere begin
using MarkovProcesses
+ import Distributed: nworkers
absolute_path = get_module_path() * "/tests/cosmos/"
# Values x1, x2 t1, t2
str_model = "ER"
@@ -32,7 +33,7 @@ for i = 1:nb_k1
k1 = l_k1[i]
k2 = l_k2[j]
command = `Cosmos $(absolute_path * "models/" * str_model * ".gspn")
- $(absolute_path * "distance_G/dist_G_" * str_model * ".lha") --njob $(ENV["JULIA_NUM_THREADS"])
+ $(absolute_path * "distance_G/dist_G_" * str_model * ".lha") --njob $(nworkers())
--const k_1=$(k1),k_2=$(k2),x1=$x1,x2=$x2,t1=$t1,t2=$t2
--level $(level) --width $(width)
--verbose 0`
@@ -47,16 +48,19 @@ for i = 1:nb_k1
set_param!(ER, "k2", convert(Float64, k2))
sync_ER = ER*A_G
mat_dist_pkg[i,j] = distribute_mean_value_lha(sync_ER, "d", nb_sim)
+ nb_accepts_pkg = distribute_prob_accept_lha(sync_ER, nb_sim)
+ #@info "About accepts" nb_sim nb_accepted nb_accepts_pkg
test = isapprox(mat_dist_cosmos[i,j], mat_dist_pkg[i,j]; atol = width*1.01)
+ test2 = nb_accepts_pkg == (nb_sim / nb_accepted)
if !test
@info "Distances too different" (k1,k2) mat_dist_pkg[i,j] mat_dist_cosmos[i,j]
end
- global test_all = test_all && test
+ global test_all = test_all && test && test2
end
end
-@info mat_dist_pkg
-@info mat_dist_cosmos
+@info "Distances R5 pkg" mat_dist_pkg
+@info "Distances R5 Cosmos" mat_dist_cosmos
rm("Result_dist_G_$(str_model).res")
rm("Result.res")
diff --git a/tests/cosmos/distance_G/dist_G_ER.lha b/tests/cosmos/distance_G/dist_G_ER.lha
index 72eeb5b..8b45dcc 100644
--- a/tests/cosmos/distance_G/dist_G_ER.lha
+++ b/tests/cosmos/distance_G/dist_G_ER.lha
@@ -47,7 +47,7 @@ Edges={
((l1,l2), #, in=1 & t>=t2, #);
((l1,l2), #, in=0 & t>=t2, {d=d*(t2-t1)});
((l1,l2), #, test_abs=0 & t<=t1, {d=d*(t2-t1)});
-((l1,l2), #, test_abs=0 & t>=t1, {d=d+(t2-t1)*min(((n-x1)^2)^0.5,((n-x2)^2)^0.5)});
+((l1,l2), #, test_abs=0 & t>=t1 & t<=t2, {d=d+(t2-t)*min(((n-x1)^2)^0.5,((n-x2)^2)^0.5)});
((l3,l1), ALL, t>=0, {n=E,test_abs=k_1*(E*S)+k_2*ES});
diff --git a/tests/run_cosmos.jl b/tests/run_cosmos.jl
index 096a7f5..bc2ad79 100644
--- a/tests/run_cosmos.jl
+++ b/tests/run_cosmos.jl
@@ -3,5 +3,6 @@ using Test
@testset "Cosmos tests" begin
@test include("cosmos/distance_F/ER_1D.jl")
+ @test include("cosmos/distance_G/ER_R5.jl")
end
--
GitLab