From c095bee4783a217aea6aa56df2e6c98e2b060178 Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Mon, 30 Nov 2020 17:21:46 +0100
Subject: [PATCH] A major fix in simulation was done, side effect issue, an
efficient copyto! is implemented for LHA states.
Automaton G works almost well, a bias in the distance is identified
compared to Cosmos when one of the model parameters is zero.
Implementation of volatile_simulate, which simulates SynchronizedModel
without saving values: it only returns the last LHA state of the
simulation.
Implementation of distributed computations of variable automata with
volatile_simulate.
Tests passed.
---
automata/automaton_F.jl | 17 +-
automata/automaton_G.jl | 158 ++++++++++++++----
automata/automaton_G_F.jl | 124 ++++++++++++++
core/MarkovProcesses.jl | 8 +-
core/lha.jl | 21 ++-
core/model.jl | 86 +++++++++-
core/plots.jl | 10 +-
core/trajectory.jl | 2 +-
tests/automata/accept_R5.jl | 29 ++++
.../automata/read_trajectory_last_state_F.jl | 6 +-
.../automata/read_trajectory_last_state_G.jl | 5 +-
tests/cosmos/distance_F/ER_1D.jl | 55 +++---
tests/cosmos/distance_G/ER_R5.jl | 61 ++++---
tests/cosmos/distance_G/dist_G_ER.lha | 19 ++-
tests/run_automata.jl | 1 +
tests/run_unit.jl | 1 +
tests/unit/models_exps_er_1d.jl | 15 ++
17 files changed, 479 insertions(+), 139 deletions(-)
create mode 100644 automata/automaton_G_F.jl
create mode 100644 tests/automata/accept_R5.jl
create mode 100644 tests/unit/models_exps_er_1d.jl
diff --git a/automata/automaton_F.jl b/automata/automaton_F.jl
index d5d452a..b464a7a 100644
--- a/automata/automaton_F.jl
+++ b/automata/automaton_F.jl
@@ -17,17 +17,19 @@ function create_automaton_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
#P <=> xn[map_var_model_idx[l_ctes[str_O]] with str_O = "P". On stock str_O dans l_ctes
# P = get_value(A, x, str_obs)
## Map of automaton variables
- map_var_automaton_idx = Dict{VariableAutomaton,Int}("n" => 1, "d" => 2)
+ map_var_automaton_idx = Dict{VariableAutomaton,Int}("n" => 1, "d" => 2, "isabs" => 3)
## 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])
+ l_flow = Dict{VariableAutomaton,Vector{Float64}}("l0" => [0.0,0.0,0.0],
+ "l1" => [0.0,0.0,0.0],
+ "l2" => [0.0,0.0,0.0],
+ "l3" => [0.0,0.0,0.0])
## Edges
map_edges = Dict{Tuple{Location,Location}, Vector{Edge}}()
+ istrue(val::Float64) = convert(Bool, val)
+
# l0 loc : we construct the edges of the form l0 => (..)
# "cc" as check_constraints
tuple = ("l0", "l1")
@@ -44,7 +46,7 @@ function create_automaton_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
us_aut_F_l1l2_1!(A::LHA, S::StateLHA, x::Vector{Int}) = (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"])
+ 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!)
@@ -76,7 +78,8 @@ function create_automaton_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
# 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::Vector{Int}) = (S["n"] = get_value(A, x, str_obs); S.loc = "l1")
+ 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))
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"])
diff --git a/automata/automaton_G.jl b/automata/automaton_G.jl
index 38f9ded..9c982c9 100644
--- a/automata/automaton_G.jl
+++ b/automata/automaton_G.jl
@@ -16,100 +16,183 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
## Map of automaton variables
map_var_automaton_idx = Dict{VariableAutomaton,Int}("t'" => 1, "in" => 2,
- "n" => 3, "d" => 4)
+ "n" => 3, "d" => 4, "isabs" => 5)
## Flow of variables
- l_flow = Dict{VariableAutomaton,Vector{Float64}}("l0" => [0.0,0.0,0.0,0.0],
- "l1" => [0.0,0.0,0.0,0.0],
- "l2" => [0.0,0.0,0.0,0.0],
- "l3" => [0.0,0.0,0.0,0.0],
- "l4" => [1.0,0.0,0.0,0.0])
+ l_flow = Dict{VariableAutomaton,Vector{Float64}}("l0" => [0.0,0.0,0.0,0.0,0.0],
+ "l1" => [0.0,0.0,0.0,0.0,0.0],
+ "l2" => [0.0,0.0,0.0,0.0,0.0],
+ "l3" => [0.0,0.0,0.0,0.0,0.0],
+ "l4" => [1.0,0.0,0.0,0.0,0.0])
## Edges
map_edges = Dict{Tuple{Location,Location}, Vector{Edge}}()
- isin(val::Float64) = convert(Bool, val)
+ istrue(val::Float64) = convert(Bool, val)
# l0 loc
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)
+ 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))
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"])); S["in"] = false)
+ (S.loc = "l3";
+ S["d"] = min(abs(A.l_ctes["x1"] - S["n"]), abs(A.l_ctes["x2"] - S["n"]));
+ S["in"] = false)
edge1 = Edge([nothing], cc_aut_G_l1l3_1, us_aut_G_l1l3_1!)
+
+ cc_aut_G_l1l3_3(A::LHA, S::StateLHA) =
+ !istrue(S["in"]) &&
+ (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"]) &&
+ (A.l_ctes["x1"] <= S["n"] <= A.l_ctes["x2"])
+ 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)
+ edge3 = Edge([nothing], cc_aut_G_l1l3_3, us_aut_G_l1l3_3!)
+
cc_aut_G_l1l3_2(A::LHA, S::StateLHA) =
- (S.time < A.l_ctes["t1"] && (A.l_ctes["x1"] <= S["n"] <= A.l_ctes["x2"]))
+ (S.time <= A.l_ctes["t1"]) &&
+ (A.l_ctes["x1"] <= S["n"] <= A.l_ctes["x2"])
us_aut_G_l1l3_2!(A::LHA, S::StateLHA, x::Vector{Int}) =
- (S.loc = "l3"; S["d"] = 0; S["in"] = false)
+ (S.loc = "l3";
+ S["d"] = 0;
+ S["in"] = false)
edge2 = Edge([nothing], cc_aut_G_l1l3_2, us_aut_G_l1l3_2!)
- cc_aut_G_l1l3_3(A::LHA, S::StateLHA) =
- (!isin(S["in"]) && (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"]) && (A.l_ctes["x1"] <= S["n"] <= A.l_ctes["x2"]))
- 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)
- edge3 = Edge([nothing], cc_aut_G_l1l3_3, us_aut_G_l1l3_3!)
+
cc_aut_G_l1l3_4(A::LHA, S::StateLHA) =
- (isin(S["in"]) && (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"]) && (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)
+ istrue(S["in"]) &&
+ (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"]) &&
+ (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)
edge4 = Edge([nothing], cc_aut_G_l1l3_4, us_aut_G_l1l3_4!)
+
map_edges[tuple] = [edge1, edge2, edge3, edge4]
tuple = ("l1", "l4")
cc_aut_G_l1l4_1(A::LHA, S::StateLHA) =
- (!isin(S["in"]) && (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"]) && (S["n"] < A.l_ctes["x1"] || S["n"] > A.l_ctes["x2"]))
- us_aut_G_l1l4_1!(A::LHA, S::StateLHA, x::Vector{Int}) = (S.loc = "l4"; S["d"] += S["d"] * (S.time - A.l_ctes["t1"]))
+ !istrue(S["in"]) &&
+ (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"]) &&
+ (S["n"] < A.l_ctes["x1"] || S["n"] > A.l_ctes["x2"])
+ us_aut_G_l1l4_1!(A::LHA, S::StateLHA, x::Vector{Int}) =
+ (S.loc = "l4";
+ S["d"] += S["d"] * (S.time - A.l_ctes["t1"]))
edge1 = Edge([nothing], cc_aut_G_l1l4_1, us_aut_G_l1l4_1!)
cc_aut_G_l1l4_2(A::LHA, S::StateLHA) =
- (isin(S["in"]) && (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"]) && (S["n"] < A.l_ctes["x1"] || S["n"] > A.l_ctes["x2"]))
- us_aut_G_l1l4_2!(A::LHA, S::StateLHA, x::Vector{Int}) = (S.loc = "l4")
+ istrue(S["in"]) &&
+ (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"]) &&
+ (S["n"] < A.l_ctes["x1"] || S["n"] > A.l_ctes["x2"])
+ us_aut_G_l1l4_2!(A::LHA, S::StateLHA, x::Vector{Int}) =
+ (S.loc = "l4")
edge2 = Edge([nothing], cc_aut_G_l1l4_2, us_aut_G_l1l4_2!)
map_edges[tuple] = [edge1, edge2]
tuple = ("l1", "l2")
- cc_aut_G_l1l2_1(A::LHA, S::StateLHA) = (isin(S["in"]) && S.time >= A.l_ctes["t2"])
- us_aut_G_l1l2_1!(A::LHA, S::StateLHA, x::Vector{Int}) = (S.loc = "l2")
+ cc_aut_G_l1l2_1(A::LHA, S::StateLHA) =
+ istrue(S["in"]) &&
+ S.time >= A.l_ctes["t2"]
+ us_aut_G_l1l2_1!(A::LHA, S::StateLHA, x::Vector{Int}) =
+ (S.loc = "l2")
edge1 = Edge([nothing], cc_aut_G_l1l2_1, us_aut_G_l1l2_1!)
- cc_aut_G_l1l2_2(A::LHA, S::StateLHA) = (!isin(S["in"]) && S.time >= A.l_ctes["t2"])
- us_aut_G_l1l2_2!(A::LHA, S::StateLHA, x::Vector{Int}) = (S.loc = "l2"; S["d"] = S["d"] * (A.l_ctes["t2"] - A.l_ctes["t1"]))
+ cc_aut_G_l1l2_2(A::LHA, S::StateLHA) =
+ !istrue(S["in"]) &&
+ S.time >= A.l_ctes["t2"]
+ us_aut_G_l1l2_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_l1l2_2, us_aut_G_l1l2_2!)
map_edges[tuple] = [edge1, edge2]
# l3 loc
tuple = ("l3", "l1")
cc_aut_G_l3l1_1(A::LHA, S::StateLHA) = true
- us_aut_G_l3l1_1!(A::LHA, S::StateLHA, x::Vector{Int}) = (S.loc = "l1"; S["n"] = get_value(A, x, str_obs))
+ us_aut_G_l3l1_1!(A::LHA, S::StateLHA, x::Vector{Int}) =
+ (S.loc = "l1";
+ S["n"] = get_value(A, x, str_obs);
+ S["isabs"] = m.isabsorbing(m.p,x))
edge1 = Edge(["ALL"], cc_aut_G_l3l1_1, us_aut_G_l3l1_1!)
map_edges[tuple] = [edge1]
tuple = ("l3", "l2")
- cc_aut_G_l3l2_1(A::LHA, S::StateLHA) = (isin(S["in"]) && S.time >= A.l_ctes["t2"])
- us_aut_G_l3l2_1!(A::LHA, S::StateLHA, x::Vector{Int}) = (S.loc = "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) = (!isin(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"]))
+ cc_aut_G_l3l2_2(A::LHA, S::StateLHA) =
+ !istrue(S["in"]) &&
+ (S.time >= A.l_ctes["t2"] || istrue(S["isabs"]))
+ 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!)
- map_edges[tuple] = [edge1, edge2]
+ 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!)
+
+ map_edges[tuple] = [edge1, edge2, edge3, edge4]
# 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["n"] = get_value(A, x, str_obs); S["in"] = true)
+ (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["n"] = get_value(A, x, str_obs);
+ S["in"] = true;
+ S["isabs"] = m.isabsorbing(m.p,x))
edge1 = Edge(["ALL"], cc_aut_G_l4l1_1, us_aut_G_l4l1_1!)
map_edges[tuple] = [edge1]
tuple = ("l4", "l2")
- cc_aut_G_l4l2_1(A::LHA, S::StateLHA) = (S.time >= A.l_ctes["t2"])
+ cc_aut_G_l4l2_1(A::LHA, S::StateLHA) =
+ (S.time >= A.l_ctes["t2"] || istrue(S["isabs"]))
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.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)
edge1 = Edge([nothing], cc_aut_G_l4l2_1, us_aut_G_l4l2_1!)
- map_edges[tuple] = [edge1]
+ 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]
## Constants
l_ctes = Dict{String,Float64}("x1" => x1, "x2" => x2, "t1" => t1, "t2" => t2)
@@ -117,6 +200,7 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
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_G
diff --git a/automata/automaton_G_F.jl b/automata/automaton_G_F.jl
new file mode 100644
index 0000000..4cd8df9
--- /dev/null
+++ b/automata/automaton_G_F.jl
@@ -0,0 +1,124 @@
+
+function create_automaton_G_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, x3::Float64, x4::Float64,
+ t1::Float64, t2::Float64, t1::Float64, t2::Float64, str_obs::String)
+ @assert str_obs in m.g
+ # Locations
+ l_loc = ["l0", "l1", "l2", "l3", "l4"]
+
+ # 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,
+ "l4" => true_inv_predicate)
+
+ ## Init and final loc
+ l_loc_init = ["l0"]
+ l_loc_final = ["l2"]
+
+ ## Map of automaton variables
+ map_var_automaton_idx = Dict{VariableAutomaton,Int}("t'" => 1, "in" => 2,
+ "n" => 3, "d" => 4)
+
+ ## Flow of variables
+ l_flow = Dict{VariableAutomaton,Vector{Float64}}("l0" => [0.0,0.0,0.0,0.0],
+ "l1" => [0.0,0.0,0.0,0.0],
+ "l2" => [0.0,0.0,0.0,0.0],
+ "l3" => [0.0,0.0,0.0,0.0],
+ "l4" => [1.0,0.0,0.0,0.0])
+
+ ## Edges
+ map_edges = Dict{Tuple{Location,Location}, Vector{Edge}}()
+
+ isin(val::Float64) = convert(Bool, val)
+
+ # l0 loc
+ 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)
+ 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"]))
+ 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"])); S["in"] = false)
+ edge1 = Edge([nothing], cc_aut_G_l1l3_1, us_aut_G_l1l3_1!)
+ cc_aut_G_l1l3_2(A::LHA, S::StateLHA) =
+ (S.time < A.l_ctes["t1"] && (A.l_ctes["x1"] <= S["n"] <= A.l_ctes["x2"]))
+ us_aut_G_l1l3_2!(A::LHA, S::StateLHA, x::Vector{Int}) =
+ (S.loc = "l3"; S["d"] = 0; S["in"] = false)
+ edge2 = Edge([nothing], cc_aut_G_l1l3_2, us_aut_G_l1l3_2!)
+ cc_aut_G_l1l3_3(A::LHA, S::StateLHA) =
+ (!isin(S["in"]) && (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"]) && (A.l_ctes["x1"] <= S["n"] <= A.l_ctes["x2"]))
+ 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)
+ edge3 = Edge([nothing], cc_aut_G_l1l3_3, us_aut_G_l1l3_3!)
+ cc_aut_G_l1l3_4(A::LHA, S::StateLHA) =
+ (isin(S["in"]) && (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"]) && (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)
+ edge4 = Edge([nothing], cc_aut_G_l1l3_4, us_aut_G_l1l3_4!)
+ map_edges[tuple] = [edge1, edge2, edge3, edge4]
+
+ tuple = ("l1", "l4")
+ cc_aut_G_l1l4_1(A::LHA, S::StateLHA) =
+ (!isin(S["in"]) && (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"]) && (S["n"] < A.l_ctes["x1"] || S["n"] > A.l_ctes["x2"]))
+ us_aut_G_l1l4_1!(A::LHA, S::StateLHA, x::Vector{Int}) = (S.loc = "l4"; S["d"] += S["d"] * (S.time - A.l_ctes["t1"]))
+ edge1 = Edge([nothing], cc_aut_G_l1l4_1, us_aut_G_l1l4_1!)
+ cc_aut_G_l1l4_2(A::LHA, S::StateLHA) =
+ (isin(S["in"]) && (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"]) && (S["n"] < A.l_ctes["x1"] || S["n"] > A.l_ctes["x2"]))
+ us_aut_G_l1l4_2!(A::LHA, S::StateLHA, x::Vector{Int}) = (S.loc = "l4")
+ edge2 = Edge([nothing], cc_aut_G_l1l4_2, us_aut_G_l1l4_2!)
+ map_edges[tuple] = [edge1, edge2]
+
+ tuple = ("l1", "l2")
+ cc_aut_G_l1l2_1(A::LHA, S::StateLHA) = (isin(S["in"]) && S.time >= A.l_ctes["t2"])
+ us_aut_G_l1l2_1!(A::LHA, S::StateLHA, x::Vector{Int}) = (S.loc = "l2")
+ edge1 = Edge([nothing], cc_aut_G_l1l2_1, us_aut_G_l1l2_1!)
+ cc_aut_G_l1l2_2(A::LHA, S::StateLHA) = (!isin(S["in"]) && S.time >= A.l_ctes["t2"])
+ us_aut_G_l1l2_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_l1l2_2, us_aut_G_l1l2_2!)
+ map_edges[tuple] = [edge1, edge2]
+
+ # l3 loc
+ tuple = ("l3", "l1")
+ cc_aut_G_l3l1_1(A::LHA, S::StateLHA) = true
+ us_aut_G_l3l1_1!(A::LHA, S::StateLHA, x::Vector{Int}) = (S.loc = "l1"; S["n"] = get_value(A, x, str_obs))
+ edge1 = Edge(["ALL"], cc_aut_G_l3l1_1, us_aut_G_l3l1_1!)
+ map_edges[tuple] = [edge1]
+
+ tuple = ("l3", "l2")
+ cc_aut_G_l3l2_1(A::LHA, S::StateLHA) = (isin(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!)
+ cc_aut_G_l3l2_2(A::LHA, S::StateLHA) = (!isin(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!)
+ 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["n"] = get_value(A, x, str_obs); S["in"] = true)
+ edge1 = Edge(["ALL"], cc_aut_G_l4l1_1, us_aut_G_l4l1_1!)
+ map_edges[tuple] = [edge1]
+
+ tuple = ("l4", "l2")
+ cc_aut_G_l4l2_1(A::LHA, S::StateLHA) = (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)
+ edge1 = Edge([nothing], cc_aut_G_l4l2_1, us_aut_G_l4l2_1!)
+ map_edges[tuple] = [edge1]
+
+ ## 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_G
+
diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index 0b4432c..ea5d129 100644
--- a/core/MarkovProcesses.jl
+++ b/core/MarkovProcesses.jl
@@ -2,10 +2,10 @@ module MarkovProcesses
import Base: +, -, *
import Base: copy, getfield, getindex, getproperty, lastindex
-import Base: setindex!, setproperty!, fill!
+import Base: setindex!, setproperty!, fill!, copyto!
import StaticArrays: SVector
-import Distributed: @everywhere
+import Distributed: @everywhere, @distributed
import Distributions: Distribution, Product, Uniform, Normal
export Distribution, Product, Uniform, Normal
@@ -26,7 +26,9 @@ export init_state, next_state!, read_trajectory
export get_index, get_value, length_var, isaccepted
# Model related methods
-export simulate, set_param!, set_time_bound!, set_observed_var!, observe_all!
+export simulate, volatile_simulate
+export distribute_mean_value_lha, mean_value_lha, distribute_prob_accept_lha
+export set_param!, set_time_bound!, set_observed_var!, observe_all!
export get_param, getproperty, get_proba_model, get_observed_var
export isbounded, isaccepted, check_consistency
export draw_model!, draw!, fill!, prior_pdf!, prior_pdf, insupport
diff --git a/core/lha.jl b/core/lha.jl
index 3d670e3..fcc3fce 100644
--- a/core/lha.jl
+++ b/core/lha.jl
@@ -15,10 +15,19 @@ function Base.show(io::IO, S::StateLHA)
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")
+ print(io, "* $var = $(S.l_var[idx]) (idx = $idx)\n")
end
end
+function Base.copyto!(Sdest::StateLHA, Ssrc::StateLHA)
+ Sdest.A = Ssrc.A
+ Sdest.loc = Ssrc.loc
+ for i = eachindex(Sdest.l_var)
+ Sdest.l_var[i] = Ssrc.l_var[i]
+ end
+ Sdest.time = Ssrc.time
+end
+
isaccepted(S::StateLHA) = (S.loc in (S.A).l_loc_final)
# Methods for synchronize / read the trajectory
@@ -72,7 +81,7 @@ 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)
+ Snplus1.l_var[i] += (A.l_flow[Sn.loc])[i]*(tnplus1 - Sn.time)
end
Snplus1.time = tnplus1
@@ -94,7 +103,10 @@ function next_state!(Snplus1::StateLHA, A::LHA,
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 break end
+ if (ind_edge == 0)
+ #if (ind_edge == 0 ||Â detected_event)
+ break
+ end
if verbose
@show first_round detected_event
@show tnplus1 tr_nplus1 xnplus1
@@ -122,6 +134,7 @@ end
# For tests purposes
function read_trajectory(A::LHA, σ::Trajectory; verbose = false)
+ @assert σ.m.d == σ.m.dobs # Model should be entirely obserbed
A_new = LHA(A, (σ.m)._map_obs_var_idx)
l_t = times(σ)
l_tr = transitions(σ)
@@ -131,10 +144,10 @@ function read_trajectory(A::LHA, σ::Trajectory; verbose = false)
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)
+ copyto!(Sn, Snplus1)
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 eb95321..da8faf3 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -46,7 +46,7 @@ function simulate(m::ContinuousTimeModel; p::Union{Nothing,AbstractVector{Float6
transitions[1] = nothing
# Values at time n
n = 1
- xn = m.x0
+ xn = copy(m.x0)
tn = m.t0
# at time n+1
isabsorbing::Bool = m.isabsorbing(p_sim,xn)
@@ -71,7 +71,7 @@ function simulate(m::ContinuousTimeModel; p::Union{Nothing,AbstractVector{Float6
break
end
n += 1
- xn = vec_x
+ copyto!(xn, vec_x)
_update_values!(full_values, times, transitions, xn, tn, l_tr[1], i)
isabsorbing = m.isabsorbing(p_sim,xn)
if isabsorbing
@@ -101,7 +101,7 @@ function simulate(m::ContinuousTimeModel; p::Union{Nothing,AbstractVector{Float6
i -= 1
break
end
- xn = vec_x
+ copyto!(xn, vec_x)
_update_values!(full_values, times, transitions,
xn, tn, l_tr[1], m.estim_min_states+size_tmp+i)
isabsorbing = m.isabsorbing(p_sim,xn)
@@ -144,7 +144,7 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
S0 = init_state(A, m.x0, m.t0)
# Values at time n
n = 1
- xn = m.x0
+ xn = copy(m.x0)
tn = m.t0
Sn = copy(S0)
isabsorbing::Bool = m.isabsorbing(p_sim,xn)
@@ -171,11 +171,11 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
break
end
n += 1
- xn = vec_x
+ copyto!(xn, vec_x)
tr_n = l_tr[1]
next_state!(Snplus1, A, xn, tn, tr_n, Sn)
_update_values!(full_values, times, transitions, xn, tn, tr_n, i)
- Sn = Snplus1
+ copyto!(Sn, Snplus1)
isabsorbing = m.isabsorbing(p_sim,xn)
isacceptedLHA = isaccepted(Snplus1)
if isabsorbing ||Â isacceptedLHA
@@ -205,12 +205,12 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
i -= 1
break
end
- xn = vec_x
+ copyto!(xn, vec_x)
tr_n = l_tr[1]
next_state!(Snplus1, A, xn, tn, tr_n, Sn)
_update_values!(full_values, times, transitions,
xn, tn, tr_n, m.estim_min_states+size_tmp+i)
- Sn = Snplus1
+ copyto!(Sn, Snplus1)
isabsorbing = m.isabsorbing(p_sim,xn)
isacceptedLHA = isaccepted(Snplus1)
if isabsorbing || isacceptedLHA
@@ -233,6 +233,54 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
return SynchronizedTrajectory(Sn, product, values, times, transitions)
end
+function volatile_simulate(product::SynchronizedModel;
+ p::Union{Nothing,AbstractVector{Float64}} = nothing, verbose::Bool = false)
+ m = product.m
+ A = product.automaton
+ p_sim = m.p
+ if p != nothing
+ p_sim = p
+ end
+ S0 = init_state(A, m.x0, m.t0)
+ # Values at time n
+ n = 1
+ xn = copy(m.x0)
+ tn = m.t0
+ 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
+ end
+ m.f!(vec_x, l_t, l_tr, xn, tn, p_sim)
+ tn = l_t[1]
+ if tn > m.time_bound
+ i -= 1
+ break
+ end
+ copyto!(xn, vec_x)
+ tr_n = l_tr[1]
+ next_state!(Snplus1, A, xn, tn, tr_n, Sn)
+ copyto!(Sn, Snplus1)
+ isabsorbing = m.isabsorbing(p_sim,xn)
+ isacceptedLHA = isaccepted(Snplus1)
+ n += 1
+ end
+ return Sn
+end
+
+
"""
`simulate(pm::ParametricModel, p_prior::AbstractVector{Float64})
@@ -246,8 +294,28 @@ function simulate(pm::ParametricModel, p_prior::AbstractVector{Float64})
return simulate(pm.m; p = full_p)
end
+"""
+ `distribute_mean_value_lha(sm::SynchronizedModel, str_var::String, nbr_stim::Int)`
+
+Distribute over workers the computation of the mean value
+of a LHA over `nbr_sim` simulations of the model.
+"""
+function distribute_mean_value_lha(sm::SynchronizedModel, str_var::String, nbr_sim::Int)
+ sum_val = @distributed (+) for i = 1:nbr_sim volatile_simulate(sm)[str_var] end
+ return sum_val / nbr_sim
+end
+
+function mean_value_lha(sm::SynchronizedModel, str_var::String, nbr_sim::Int)
+ sum_val = 0.0
+ for i = 1:nbr_sim
+ sum_val += volatile_simulate(sm)[str_var]
+ end
+ return sum_val / nbr_sim
+end
-function simulate(m::ContinuousTimeModel, n::Int)
+function distribute_prob_accept_lha(sm::SynchronizedModel, str_var::String, nbr_sim::Int)
+ sum_val = @distributed (+) for i = 1:nbr_sim Int(isaccepted(volatile_simulate(sm))) end
+ return sum_val / nbr_sim
end
isbounded(m::Model) = get_proba_model(m).time_bound < Inf
diff --git a/core/plots.jl b/core/plots.jl
index d174d15..9e2ac46 100644
--- a/core/plots.jl
+++ b/core/plots.jl
@@ -1,5 +1,5 @@
-import Plots: plot, plot!, scatter!
+import Plots: plot, plot!, scatter!, Shape
import Plots: current, palette, display, png, close
"""
@@ -50,5 +50,11 @@ function plot(σ::AbstractTrajectory, vars::String...; filename::String = "", pl
end
end
-export plot
+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)
+end
+
+export plot, plot!
diff --git a/core/trajectory.jl b/core/trajectory.jl
index 730340c..89a7971 100644
--- a/core/trajectory.jl
+++ b/core/trajectory.jl
@@ -162,7 +162,6 @@ function check_consistency(σ::AbstractTrajectory)
@assert length_obs_var(σ) == get_proba_model(σ.m).dobs
return true
end
-issteadystate(σ::AbstractTrajectory) = get_proba_model(σ.m).isabsorbing(get_proba_model(σ.m).p, view(reshape(σ[end], 1, m.d), 1, :))
# Properties of the trajectory
length_states(σ::AbstractTrajectory) = length(σ.times)
@@ -170,6 +169,7 @@ length_obs_var(σ::AbstractTrajectory) = length(σ.values)
get_obs_var(σ::AbstractTrajectory) = get_proba_model(σ.m).g
isbounded(σ::AbstractTrajectory) = σ.transitions[end] == nothing && length_states(σ) >= 2
isaccepted(σ::SynchronizedTrajectory) = isaccepted(σ.S)
+issteadystate(σ::AbstractTrajectory) = @warn "Unimplemented"
# Access to trajectory values
get_var_values(σ::AbstractTrajectory, var::String) = σ.values[get_proba_model(σ.m)._map_obs_var_idx[var]]
diff --git a/tests/automata/accept_R5.jl b/tests/automata/accept_R5.jl
new file mode 100644
index 0000000..a9ab3fb
--- /dev/null
+++ b/tests/automata/accept_R5.jl
@@ -0,0 +1,29 @@
+
+using MarkovProcesses
+# Values x1, x2 t1, t2
+str_model = "ER"
+load_model(str_model)
+model = ER
+observe_all!(ER)
+ER.buffer_size = 100
+load_automaton("automaton_G")
+width = 0.5
+level = 0.95
+x1, x2, t1, t2 = 50.0, 100.0, 0.0, 0.8
+A_G = create_automaton_G(model, x1, x2, t1, t2, "P")
+set_param!(ER, "k1", 0.2)
+set_param!(ER, "k2", 40.0)
+
+test_all = true
+sync_ER = ER*A_G
+σ = nothing
+nb_sim = 1000
+for i = 1:nb_sim
+ global σ = simulate(sync_ER)
+ if !isaccepted(σ)
+ error("Ouch")
+ end
+end
+
+return test_all
+
diff --git a/tests/automata/read_trajectory_last_state_F.jl b/tests/automata/read_trajectory_last_state_F.jl
index 17f02b8..c1c799f 100644
--- a/tests/automata/read_trajectory_last_state_F.jl
+++ b/tests/automata/read_trajectory_last_state_F.jl
@@ -4,6 +4,7 @@ using MarkovProcesses
load_model("SIR")
load_automaton("automaton_F")
SIR.time_bound = 120.0
+observe_all!(SIR)
x1, x2, t1, t2 = 0.0, Inf, 100.0, 120.0
A_F = create_automaton_F(SIR, x1, x2, t1, t2, "I") # <: LHA
@@ -11,13 +12,12 @@ 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)
- #=
+ test = (get_state_from_time(σ, Send.time)[2] == Send["n"]) && (Send["d"] == 0)
if !test
@show Send
@show get_state_from_time(σ, Send.time)
error("tkt")
- end=#
+ end
return test
end
diff --git a/tests/automata/read_trajectory_last_state_G.jl b/tests/automata/read_trajectory_last_state_G.jl
index 305448b..cdd5a43 100644
--- a/tests/automata/read_trajectory_last_state_G.jl
+++ b/tests/automata/read_trajectory_last_state_G.jl
@@ -3,6 +3,7 @@ using MarkovProcesses
load_model("ER")
load_automaton("automaton_G")
+observe_all!(ER)
ER.time_bound = 2.0
x1, x2, t1, t2 = 0.0, Inf, 0.0, 2.0
@@ -11,14 +12,12 @@ A_G = create_automaton_G(ER, x1, x2, t1, t2, "P") # <: 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)
- #=
+ test = (get_state_from_time(σ, Send.time)[4] == Send["n"]) && (Send["d"] == 0)
if !test
@show Send
@show get_state_from_time(σ, Send.time)
error("bad")
end
- =#
return test
end
diff --git a/tests/cosmos/distance_F/ER_1D.jl b/tests/cosmos/distance_F/ER_1D.jl
index 0b2cc8c..80c8796 100644
--- a/tests/cosmos/distance_F/ER_1D.jl
+++ b/tests/cosmos/distance_F/ER_1D.jl
@@ -1,28 +1,26 @@
-
-using MarkovProcesses
-absolute_path = get_module_path() * "/tests/cosmos/"
-
-# Remarque: la valeur estimée par Cosmos varie plus que de 0.01
-
-# Values x1, x2 t1, t2
-dict_exp = Dict(
- "R1" => [50, 75, 0.025, 0.05],
- "R2" => [50, 75, 0.05, 0.075],
- "R3" => [25, 50, 0.05, 0.075]
- )
-
-str_model = "ER"
-if str_model == "ER"
- load_model(str_model)
- model = ER
+@everywhere begin
+ using MarkovProcesses
+ import Distributed: nworkers
+ absolute_path = get_module_path() * "/tests/cosmos/"
+ # Values x1, x2 t1, t2
+ dict_exp = Dict(
+ "R1" => [50, 75, 0.025, 0.05],
+ "R2" => [50, 75, 0.05, 0.075],
+ "R3" => [25, 50, 0.05, 0.075]
+ )
+ str_model = "ER"
+ if str_model == "ER"
+ load_model(str_model)
+ model = ER
+ observe_all!(ER)
+ end
+ load_automaton("automaton_F")
+ test_all = true
+ width = 0.2
+ level = 0.95
+ l_exp = ["R1", "R2", "R3"]
end
-load_automaton("automaton_F")
-
-test_all = true
-width = 0.1
-level = 0.95
-l_exp = ["R1", "R2", "R3"]
#l_exp = ["R2"]
for exp in l_exp
if !(exp in keys(dict_exp))
@@ -40,7 +38,7 @@ for exp in l_exp
# Cosmos estimation
k3 = l_k3[i]
command = `Cosmos $(absolute_path * "models/" * str_model * ".gspn")
- $(absolute_path * "distance_F/dist_F_" * str_model * ".lha") --njob $(ENV["JULIA_NUM_THREADS"])
+ $(absolute_path * "distance_F/dist_F_" * str_model * ".lha") --njob $(nworkers())
--const k_3=$(k3),x1=$x1,x2=$x2,t1=$t1,t2=$t2
--level $(level) --width $(width)
--verbose 0`
@@ -53,20 +51,17 @@ for exp in l_exp
#Â MarkovProcesses estimation
set_param!(ER, "k3", convert(Float64, k3))
sync_ER = ER*A_F
- l_dist_i_threads = zeros(Threads.nthreads())
- Threads.@threads for j = 1:nb_sim
- σ = simulate(sync_ER)
- l_dist_i_threads[Threads.threadid()] += (σ.S)["d"]
- end
- l_dist_pkg[i] = sum(l_dist_i_threads) / nb_sim
+ l_dist_pkg[i] = distribute_mean_value_lha(sync_ER, "d", nb_sim)
test = isapprox(l_dist_cosmos[i], l_dist_pkg[i]; atol = width*1.01)
global test_all = test_all && test
if !test
@show l_dist_pkg[i], l_dist_cosmos[i]
@show exp
@show k3
+ @show nb_sim
end
end
+ @info exp l_dist_pkg l_dist_cosmos
end
rm("Result_dist_F_$(str_model).res")
diff --git a/tests/cosmos/distance_G/ER_R5.jl b/tests/cosmos/distance_G/ER_R5.jl
index c9655ec..b211a29 100644
--- a/tests/cosmos/distance_G/ER_R5.jl
+++ b/tests/cosmos/distance_G/ER_R5.jl
@@ -1,26 +1,25 @@
-using MarkovProcesses
-absolute_path = get_module_path() * "/tests/cosmos/"
-
-# Remarque: la valeur estimée par Cosmos varie plus que de 0.01
-
-# Values x1, x2 t1, t2
-str_model = "ER"
-load_model(str_model)
-model = ER
-ER.buffer_size = 100
-load_automaton("automaton_G")
+@everywhere begin
+ using MarkovProcesses
+ absolute_path = get_module_path() * "/tests/cosmos/"
+ # Values x1, x2 t1, t2
+ str_model = "ER"
+ load_model(str_model)
+ model = ER
+ observe_all!(ER)
+ ER.buffer_size = 100
+ load_automaton("automaton_G")
+ width = 0.5
+ level = 0.95
+ x1, x2, t1, t2 = 50.0, 100.0, 0.0, 0.8
+ A_G = create_automaton_G(model, x1, x2, t1, t2, "E")
+ l_k1 = 0.0:0.5:1.5
+ #l_k1 = 0.2:0.2
+ l_k2 = 0:40:100
+ #l_k2 = 40:40
+end
test_all = true
-width = 0.1
-level = 0.95
-x1, x2, t1, t2 = 50.0, 100.0, 0.0, 0.8
-A_G = create_automaton_G(model, x1, x2, t1, t2, "P")
-
-l_k1 = 0.0:0.2:1.4
-l_k1 = 0.2:0.2
-l_k2 = 0:20:100
-l_k2 = 40:40
nb_k1 = length(l_k1)
nb_k2 = length(l_k2)
mat_dist_cosmos = zeros(nb_k1,nb_k2)
@@ -34,11 +33,10 @@ for i = 1:nb_k1
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"])
- --const k_1=$(k1),k2=$(k2),x1=$x1,x2=$x2,t1=$t1,t2=$t2
+ --const k_1=$(k1),k_2=$(k2),x1=$x1,x2=$x2,t1=$t1,t2=$t2
--level $(level) --width $(width)
- --verbose 2`
- #run(pipeline(command, stderr=devnull))
- @timev run(pipeline(command))
+ --verbose 0`
+ run(pipeline(command, stderr=devnull))
dict_values = cosmos_get_values("Result_dist_G_$(str_model).res")
mat_dist_cosmos[i,j] = dict_values["Estimated value"]
nb_sim = dict_values["Total paths"]
@@ -48,18 +46,17 @@ for i = 1:nb_k1
set_param!(ER, "k1", convert(Float64, k1))
set_param!(ER, "k2", convert(Float64, k2))
sync_ER = ER*A_G
- mat_dist_ij_threads = zeros(Threads.nthreads())
- @timev Threads.@threads for j = 1:nb_sim
- σ = simulate(sync_ER)
- mat_dist_ij_threads[Threads.threadid()] += (σ.S)["d"]
+ mat_dist_pkg[i,j] = distribute_mean_value_lha(sync_ER, "d", nb_sim)
+ test = isapprox(mat_dist_cosmos[i,j], mat_dist_pkg[i,j]; atol = width*1.01)
+ if !test
+ @info "Distances too different" (k1,k2) mat_dist_pkg[i,j] mat_dist_cosmos[i,j]
end
- mat_dist_pkg[i,j] = sum(mat_dist_ij_threads) / nb_sim
- global test_all = test_all && isapprox(mat_dist_cosmos[i,j], mat_dist_pkg[i,j]; atol = width)
+ global test_all = test_all && test
end
end
-@show mat_dist_pkg
-@show mat_dist_cosmos
+@info mat_dist_pkg
+@info 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 13872bf..72eeb5b 100644
--- a/tests/cosmos/distance_G/dist_G_ER.lha
+++ b/tests/cosmos/distance_G/dist_G_ER.lha
@@ -31,28 +31,31 @@ Locations={
Edges={
% Init
-((l0,l1), #, t>=0, {n=E,d=0,in=1,test_abs=(E*S)+ES});
+((l0,l1), #, t>=0, {n=E,d=0,in=1,test_abs=k_1*(E*S)+k_2*ES});
((l1,l3), #, t<=t1 & n<=x1-1, {d=min(((n-x1)^2)^0.5,((n-x2)^2)^0.5),in=0});
((l1,l3), #, t<=t1 & n>=x2+1, {d=min(((n-x1)^2)^0.5,((n-x2)^2)^0.5),in=0});
((l1,l3), #, in=0 & t>=t1 & t<=t2 & n>=x1 & n<=x2, {d=d*(t-t1),tprime=0});
((l1,l3), #, t<=t1 & n>=x1 & n<=x2, {d=0,in=0});
-((l1,l3), #, in=1 & t>=t1 & t<=t2 & n>=x1 & n<=x2, {tprime=0,in=0});
-((l3,l1), ALL, t>=0, {n=E,test_abs=(E*S)+ES});
+((l1,l3), #, in=1 & t>=t1 & t<=t2 & n>=x1 & n<=x2, {tprime=0});
((l1,l4), #, in=0 & t>=t1 & t<=t2 & n<=x1-1, {d=d+d*(t-t1)});
((l1,l4), #, in=0 & t>=t1 & t<=t2 & n>=x2+1, {d=d+d*(t-t1)});
((l1,l4), #, in=1 & t>=t1 & t<=t2 & n<=x1-1, #);
((l1,l4), #, in=1 & t>=t1 & t<=t2 & n>=x2+1, #);
-((l4,l1), ALL, t>=0, {d=d+tprime*min(((n-x1)^2)^0.5,((n-x2)^2)^0.5),tprime=0,n=E,in=1,test_abs=(E*S)+ES});
-% Acceptation
-((l4,l2), #, t>=t2, {d=d+tprime*min(((n-x1)^2)^0.5,((n-x2)^2)^0.5)});
-((l3,l2), #, in=1 & t>=t2, {d=d*(t2-t1)});
-((l3,l2), #, in=0 & t>=t2, #);
((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)});
+
+((l3,l1), ALL, t>=0, {n=E,test_abs=k_1*(E*S)+k_2*ES});
+
+((l3,l2), #, in=1 & t>=t2, {d=d*(t2-t1)});
+((l3,l2), #, in=0 & t>=t2, #);
+
+((l4,l1), ALL, t>=0, {d=d+tprime*min(((n-x1)^2)^0.5,((n-x2)^2)^0.5),tprime=0,n=E,in=1,test_abs=k_1*(E*S)+k_2*ES});
+
+((l4,l2), #, t>=t2, {d=d+tprime*min(((n-x1)^2)^0.5,((n-x2)^2)^0.5)});
};
diff --git a/tests/run_automata.jl b/tests/run_automata.jl
index ef0e2ee..e94f7c5 100644
--- a/tests/run_automata.jl
+++ b/tests/run_automata.jl
@@ -6,5 +6,6 @@ using Test
@test include("automata/read_trajectory_last_state_G.jl")
@test include("automata/sync_simulate_last_state_F.jl")
@test include("automata/sync_simulate_last_state_G.jl")
+ @test include("automata/accept_R5.jl")
end
diff --git a/tests/run_unit.jl b/tests/run_unit.jl
index 46dc9d9..1fbcbf6 100644
--- a/tests/run_unit.jl
+++ b/tests/run_unit.jl
@@ -27,6 +27,7 @@ using Test
@test include("unit/long_sim_er.jl")
@test include("unit/model_prior.jl")
+ @test include("unit/models_exps_er_1d.jl")
@test include("unit/observe_all.jl")
diff --git a/tests/unit/models_exps_er_1d.jl b/tests/unit/models_exps_er_1d.jl
new file mode 100644
index 0000000..704516b
--- /dev/null
+++ b/tests/unit/models_exps_er_1d.jl
@@ -0,0 +1,15 @@
+
+using MarkovProcesses
+
+load_model("ER")
+observe_all!(ER)
+load_automaton("automaton_F")
+load_automaton("automaton_G")
+
+A_F_R1 = create_automaton_F(ER, 50.0, 75.0, 0.025, 0.05, "P")
+A_F_R2 = create_automaton_F(ER, 50.0, 75.0, 0.05, 0.075, "P")
+A_F_R3 = create_automaton_F(ER, 25.0, 50.0, 0.05, 0.075, "P")
+A_G_R5 = create_automaton_G(ER, 50.0, 100.0, 0.0, 0.8, "E")
+
+return true
+
--
GitLab