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,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"])

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

automata/automaton_G_F.jl

+
+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
+

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

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

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