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

@@ -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]
 

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)

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
 

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

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!

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

tests/cosmos/distance_F/dist_F_ER.lha

-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, #);
 };
+

tests/cosmos/distance_G/ER_R5.jl

 
 @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")

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});
 

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