Skip to content
Snippets Groups Projects
Commit 146e2afa authored by Bentriou Mahmoud's avatar Bentriou Mahmoud
Browse files

Improvement of the simulation core, less memory allocated.

parent 618d7910
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 177 additions and 248 deletions
automata/automaton_F.jl
+ 9
− 9
View file @ 146e2afa
...@@ -32,7 +32,7 @@ function create_automaton_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1 ...@@ -32,7 +32,7 @@ function create_automaton_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
# "cc" as check_constraints # "cc" as check_constraints
tuple = ("l0", "l1") tuple = ("l0", "l1")
cc_aut_F_l0l1_1(A::LHA, S::StateLHA) = true cc_aut_F_l0l1_1(A::LHA, S::StateLHA) = true
us_aut_F_l0l1_1!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = 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))
edge1 = Edge([nothing], cc_aut_F_l0l1_1, us_aut_F_l0l1_1!) edge1 = Edge([nothing], cc_aut_F_l0l1_1, us_aut_F_l0l1_1!)
map_edges[tuple] = [edge1] map_edges[tuple] = [edge1]
...@@ -41,34 +41,34 @@ function create_automaton_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1 ...@@ -41,34 +41,34 @@ function create_automaton_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
tuple = ("l1", "l2") tuple = ("l1", "l2")
cc_aut_F_l1l2_1(A::LHA, S::StateLHA) = 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"]) (A.l_ctes["x1"] <= S["n"] <= A.l_ctes["x2"]) && (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"])
us_aut_F_l1l2_1!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = (S["d"] = 0; S.loc = "l2") 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!) 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"])
us_aut_F_l1l2_2!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = (S.loc = "l2") 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!) 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"]) cc_aut_F_l1l2_3(A::LHA, S::StateLHA) = (S["d"] == 0 && S.time >= A.l_ctes["t1"])
us_aut_F_l1l2_3!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = (S.loc = "l2") 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!) edge3 = Edge([nothing], cc_aut_F_l1l2_3, us_aut_F_l1l2_3!)
map_edges[tuple] = [edge1, edge2, edge3] map_edges[tuple] = [edge1, edge2, edge3]
tuple = ("l1", "l3") tuple = ("l1", "l3")
cc_aut_F_l1l3_1(A::LHA, S::StateLHA) = (A.l_ctes["x1"] <= S["n"] <= A.l_ctes["x2"]) cc_aut_F_l1l3_1(A::LHA, S::StateLHA) = (A.l_ctes["x1"] <= S["n"] <= A.l_ctes["x2"])
us_aut_F_l1l3_1!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = (S["d"] = 0; S.loc = "l3") us_aut_F_l1l3_1!(A::LHA, S::StateLHA, x::Vector{Int}) = (S["d"] = 0; S.loc = "l3")
edge1 = Edge([nothing], cc_aut_F_l1l3_1, us_aut_F_l1l3_1!) edge1 = Edge([nothing], cc_aut_F_l1l3_1, us_aut_F_l1l3_1!)
cc_aut_F_l1l3_2(A::LHA, S::StateLHA) = 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"]) (A.l_ctes["x1"] > S["n"] || S["n"] > A.l_ctes["x2"]) && (S.time < A.l_ctes["t1"])
us_aut_F_l1l3_2!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = 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), (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") sqrt((S.time - A.l_ctes["t1"])^2 + (S["n"] - A.l_ctes["x1"])^2)); S.loc = "l3")
edge2 = Edge([nothing], cc_aut_F_l1l3_2, us_aut_F_l1l3_2!) edge2 = Edge([nothing], cc_aut_F_l1l3_2, us_aut_F_l1l3_2!)
cc_aut_F_l1l3_3(A::LHA, S::StateLHA) = 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"]) (A.l_ctes["x1"] > S["n"] || S["n"] > A.l_ctes["x2"]) && (A.l_ctes["t1"] <= S.time <= A.l_ctes["t2"])
us_aut_F_l1l3_3!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = 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["d"] = min(S["d"], min(abs(S["n"] - A.l_ctes["x1"]), abs(S["n"] - A.l_ctes["x2"]))); S.loc = "l3")
edge3 = Edge([nothing], cc_aut_F_l1l3_3, us_aut_F_l1l3_3!) edge3 = Edge([nothing], cc_aut_F_l1l3_3, us_aut_F_l1l3_3!)
map_edges[tuple] = [edge1, edge2, edge3] map_edges[tuple] = [edge1, edge2, edge3]
...@@ -76,11 +76,11 @@ function create_automaton_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1 ...@@ -76,11 +76,11 @@ function create_automaton_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
# l3 loc # l3 loc
tuple = ("l3", "l1") tuple = ("l3", "l1")
cc_aut_F_l3l1_1(A::LHA, S::StateLHA) = true cc_aut_F_l3l1_1(A::LHA, S::StateLHA) = true
us_aut_F_l3l1_1!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = (S["n"] = get_value(A, x, str_obs); S.loc = "l1") us_aut_F_l3l1_1!(A::LHA, S::StateLHA, x::Vector{Int}) = (S["n"] = get_value(A, x, str_obs); S.loc = "l1")
edge1 = Edge(["ALL"], cc_aut_F_l3l1_1, us_aut_F_l3l1_1!) edge1 = Edge(["ALL"], cc_aut_F_l3l1_1, us_aut_F_l3l1_1!)
tuple = ("l3", "l2") tuple = ("l3", "l2")
cc_aut_F_l3l2_1(A::LHA, S::StateLHA) = (S.time >= A.l_ctes["t2"]) cc_aut_F_l3l2_1(A::LHA, S::StateLHA) = (S.time >= A.l_ctes["t2"])
us_aut_F_l3l2_1!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = (S["n"] = get_value(A, x, str_obs); S.loc = "l2") us_aut_F_l3l2_1!(A::LHA, S::StateLHA, x::Vector{Int}) = (S["n"] = get_value(A, x, str_obs); S.loc = "l2")
edge2 = Edge([nothing], cc_aut_F_l3l2_1, us_aut_F_l3l2_1!) edge2 = Edge([nothing], cc_aut_F_l3l2_1, us_aut_F_l3l2_1!)
map_edges[tuple] = [edge1, edge2] map_edges[tuple] = [edge1, edge2]
......
automata/automaton_G.jl
+ 14
− 14
View file @ 146e2afa
...@@ -33,7 +33,7 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1 ...@@ -33,7 +33,7 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
# l0 loc # l0 loc
tuple = ("l0", "l1") tuple = ("l0", "l1")
cc_aut_G_l0l1_1(A::LHA, S::StateLHA) = true cc_aut_G_l0l1_1(A::LHA, S::StateLHA) = true
us_aut_G_l0l1_1!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = (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)
edge1 = Edge([nothing], cc_aut_G_l0l1_1, us_aut_G_l0l1_1!) edge1 = Edge([nothing], cc_aut_G_l0l1_1, us_aut_G_l0l1_1!)
map_edges[tuple] = [edge1] map_edges[tuple] = [edge1]
...@@ -41,64 +41,64 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1 ...@@ -41,64 +41,64 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
tuple = ("l1", "l3") tuple = ("l1", "l3")
cc_aut_G_l1l3_1(A::LHA, S::StateLHA) = 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::SubArray{Int,1}) = 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!) edge1 = Edge([nothing], cc_aut_G_l1l3_1, us_aut_G_l1l3_1!)
cc_aut_G_l1l3_2(A::LHA, S::StateLHA) = 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::SubArray{Int,1}) = 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!) edge2 = Edge([nothing], cc_aut_G_l1l3_2, us_aut_G_l1l3_2!)
cc_aut_G_l1l3_3(A::LHA, S::StateLHA) = 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"])) (!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::SubArray{Int,1}) = (S.loc = "l3"; S["d"] = S["d"] * (S.time - A.l_ctes["t1"]); S["t'"] = 0.0) 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!) edge3 = Edge([nothing], cc_aut_G_l1l3_3, us_aut_G_l1l3_3!)
cc_aut_G_l1l3_4(A::LHA, S::StateLHA) = 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"])) (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::SubArray{Int,1}) = (S.loc = "l3"; S["t'"] = 0.0) 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!) edge4 = Edge([nothing], cc_aut_G_l1l3_4, us_aut_G_l1l3_4!)
map_edges[tuple] = [edge1, edge2, edge3, edge4] map_edges[tuple] = [edge1, edge2, edge3, edge4]
tuple = ("l1", "l4") tuple = ("l1", "l4")
cc_aut_G_l1l4_1(A::LHA, S::StateLHA) = 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"])) (!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::SubArray{Int,1}) = (S.loc = "l4"; S["d"] += S["d"] * (S.time - A.l_ctes["t1"])) 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!) edge1 = Edge([nothing], cc_aut_G_l1l4_1, us_aut_G_l1l4_1!)
cc_aut_G_l1l4_2(A::LHA, S::StateLHA) = 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"])) (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::SubArray{Int,1}) = (S.loc = "l4") 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!) edge2 = Edge([nothing], cc_aut_G_l1l4_2, us_aut_G_l1l4_2!)
map_edges[tuple] = [edge1, edge2] map_edges[tuple] = [edge1, edge2]
tuple = ("l1", "l2") tuple = ("l1", "l2")
cc_aut_G_l1l2_1(A::LHA, S::StateLHA) = (isin(S["in"]) && S.time >= A.l_ctes["t2"]) 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::SubArray{Int,1}) = (S.loc = "l2") 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!) 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"]) 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::SubArray{Int,1}) = (S.loc = "l2"; S["d"] = S["d"] * (A.l_ctes["t2"] - A.l_ctes["t1"])) 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!) edge2 = Edge([nothing], cc_aut_G_l1l2_2, us_aut_G_l1l2_2!)
map_edges[tuple] = [edge1, edge2] map_edges[tuple] = [edge1, edge2]
# l3 loc # l3 loc
tuple = ("l3", "l1") tuple = ("l3", "l1")
cc_aut_G_l3l1_1(A::LHA, S::StateLHA) = true cc_aut_G_l3l1_1(A::LHA, S::StateLHA) = true
us_aut_G_l3l1_1!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = (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))
edge1 = Edge(["ALL"], cc_aut_G_l3l1_1, us_aut_G_l3l1_1!) edge1 = Edge(["ALL"], cc_aut_G_l3l1_1, us_aut_G_l3l1_1!)
map_edges[tuple] = [edge1] map_edges[tuple] = [edge1]
tuple = ("l3", "l2") tuple = ("l3", "l2")
cc_aut_G_l3l2_1(A::LHA, S::StateLHA) = (isin(S["in"]) && S.time >= A.l_ctes["t2"]) 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::SubArray{Int,1}) = (S.loc = "l2") 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!) 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"]) 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::SubArray{Int,1}) = (S.loc = "l2"; S["d"] = S["d"] * (A.l_ctes["t2"] - A.l_ctes["t1"])) 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!) edge2 = Edge([nothing], cc_aut_G_l3l2_2, us_aut_G_l3l2_2!)
map_edges[tuple] = [edge1, edge2] map_edges[tuple] = [edge1, edge2]
# l4 loc # l4 loc
tuple = ("l4", "l1") tuple = ("l4", "l1")
cc_aut_G_l4l1_1(A::LHA, S::StateLHA) = true cc_aut_G_l4l1_1(A::LHA, S::StateLHA) = true
us_aut_G_l4l1_1!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = 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.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["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!) edge1 = Edge(["ALL"], cc_aut_G_l4l1_1, us_aut_G_l4l1_1!)
...@@ -106,7 +106,7 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1 ...@@ -106,7 +106,7 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
tuple = ("l4", "l2") 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"])
us_aut_G_l4l2_1!(A::LHA, S::StateLHA, x::SubArray{Int,1}) = 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!) edge1 = Edge([nothing], cc_aut_G_l4l2_1, us_aut_G_l4l2_1!)
map_edges[tuple] = [edge1] map_edges[tuple] = [edge1]
......
bench/pygmalion/multiple_long_sim_er.jl
+ 8
− 8
View file @ 146e2afa
...@@ -20,29 +20,29 @@ g_all = create_observation_function([ObserverModel(str_oml, tml)]) ...@@ -20,29 +20,29 @@ g_all = create_observation_function([ObserverModel(str_oml, tml)])
@timev pygmalion.simulate(f, g_all, x0, u, p_true; on = nothing, full_timeline = true) @timev pygmalion.simulate(f, g_all, x0, u, p_true; on = nothing, full_timeline = true)
b1_pyg = @benchmark for i=1:$(nb_sim) pygmalion.simulate($f, $g_all, $x0, $u, $p_true; on = nothing, full_timeline = true) end b1_pyg = @benchmark for i=1:$(nb_sim) pygmalion.simulate($f, $g_all, $x0, $u, $p_true; on = nothing, full_timeline = true) end
b2_pyg = @benchmark for i=1:$(nb_sim) pygmalion.simulate(f, g_all, x0, u, p_true; on = nothing, full_timeline = true) end #b2_pyg = @benchmark for i=1:$(nb_sim) pygmalion.simulate(f, g_all, x0, u, p_true; on = nothing, full_timeline = true) end
@show minimum(b1_pyg), mean(b1_pyg), maximum(b1_pyg) @show minimum(b1_pyg), mean(b1_pyg), maximum(b1_pyg)
@show minimum(b2_pyg), mean(b2_pyg), maximum(b2_pyg) #@show minimum(b2_pyg), mean(b2_pyg), maximum(b2_pyg)
println("MarkovProcesses:") println("MarkovProcesses:")
using MarkovProcesses using MarkovProcesses
MarkovProcesses.load_model("ER") MarkovProcesses.load_model("ER")
ER.time_bound = 20.0 ER.time_bound = 20.0
ER.estim_min_states = 9000
set_param!(ER, "k1", 0.2) set_param!(ER, "k1", 0.2)
set_param!(ER, "k2", 40.0) set_param!(ER, "k2", 40.0)
@timev MarkovProcesses.simulate(ER) @timev MarkovProcesses.simulate(ER)
println("Default buffer size=10") println("Default buffer size=10")
b1 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate($ER) end b1 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate($ER) end
b2 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end #b2 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end
@show minimum(b1), mean(b1), maximum(b1) @show minimum(b1), mean(b1), maximum(b1)
@show minimum(b2), mean(b2), maximum(b2) #@show minimum(b2), mean(b2), maximum(b2)
println("Buffer size 100") println("Buffer size 100 / min_states 8000")
ER.estim_min_states = 8000
ER.buffer_size = 100 ER.buffer_size = 100
b1_100 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate($ER) end b1_100 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate($ER) end
b2_100 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end #b2_100 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end
@show minimum(b1_100), mean(b1_100), maximum(b1_100) @show minimum(b1_100), mean(b1_100), maximum(b1_100)
@show minimum(b2_100), mean(b2_100), maximum(b2_100) #@show minimum(b2_100), mean(b2_100), maximum(b2_100)
core/_tests_simulate.jl
+ 0
− 63
View file @ 146e2afa
...@@ -334,66 +334,3 @@ function _simulate_d7458(m::ContinuousTimeModel) ...@@ -334,66 +334,3 @@ function _simulate_d7458(m::ContinuousTimeModel)
return Trajectory(m, values, times, transitions) return Trajectory(m, values, times, transitions)
end end
function _simulate_last(product::SynchronizedModel)
# trajectory fields
m = product.m
A = product.automaton
full_values = Vector{Vector{Int}}(undef, m.d)
for i = 1:m.d full_values[i] = Int[m.x0[i]] end
for i = 1:m.d sizehint!(full_values[i], m.estim_min_states) end
times = Float64[m.t0]
transitions = Union{String,Nothing}[nothing]
reshaped_x0 = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
S0 = init_state(A, m.x0, m.t0)
# values at time n
n = 0
xn = reshaped_x0
tn = m.t0
Sn = copy(S0)
# at time n+1
mat_x = zeros(Int, m.buffer_size, m.d)
l_t = zeros(Float64, m.buffer_size)
l_tr = Vector{String}(undef, m.buffer_size)
isabsorbing::Bool = m.isabsorbing(m.p,xn)
isacceptedLHA::Bool = isaccepted(Sn)
Snplus1 = copy(Sn)
while !isabsorbing && tn <= m.time_bound && !isacceptedLHA
i = 0
while i < m.buffer_size && !isabsorbing && tn <= m.time_bound && !isacceptedLHA
i += 1
m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
tn = l_t[i]
if tn > m.time_bound
i -= 1 # 0 is an ok value, 1:0 is allowed
break
end
xn = view(mat_x, i, :)
tr_n = l_tr[i]
next_state!(Snplus1, A, xn, tn, tr_n, Sn)
Sn = Snplus1
isabsorbing = m.isabsorbing(m.p,xn)
isacceptedLHA = isaccepted(Snplus1)
end
for k = 1:m.d
append!(full_values[k], view(mat_x, 1:i, k))
end
append!(times, view(l_t, 1:i))
append!(transitions, view(l_tr, 1:i))
n += i
end
# When the trajectory is accepted, we should not add an end value
if isbounded(m) && !isaccepted(Sn)
@assert times[end] < m.time_bound
for k = 1:m.d
push!(full_values[k], full_values[k][end])
end
push!(times, m.time_bound)
push!(transitions, nothing)
end
values = full_values[m._g_idx]
return SynchronizedTrajectory(Snplus1, product, values, times, transitions)
end
core/lha.jl
+ 3
− 8
View file @ 146e2afa
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
load_automaton(automaton::String) = include(get_module_path() * "/automata/$(automaton).jl") load_automaton(automaton::String) = include(get_module_path() * "/automata/$(automaton).jl")
length_var(A::LHA) = length(A.map_var_automaton_idx) length_var(A::LHA) = length(A.map_var_automaton_idx)
get_value(A::LHA, x::SubArray{Int,1}, var::String) = x[A.map_var_model_idx[var]] get_value(A::LHA, x::Vector{Int}, var::String) = x[A.map_var_model_idx[var]]
copy(S::StateLHA) = StateLHA(S.A, S.loc, S.l_var, S.time) copy(S::StateLHA) = StateLHA(S.A, S.loc, S.l_var, S.time)
# Not overring getproperty, setproperty to avoid a conversion Symbol => String for the dict key # Not overring getproperty, setproperty to avoid a conversion Symbol => String for the dict key
...@@ -36,7 +36,7 @@ function init_state(A::LHA, x0::Vector{Int}, t0::Float64) ...@@ -36,7 +36,7 @@ function init_state(A::LHA, x0::Vector{Int}, t0::Float64)
end end
function next_state!(Snplus1::StateLHA, A::LHA, function next_state!(Snplus1::StateLHA, A::LHA,
xnplus1::SubArray{Int,1}, tnplus1::Float64, tr_nplus1::Transition, xnplus1::Vector{Int}, tnplus1::Float64, tr_nplus1::Transition,
Sn::StateLHA; verbose = false) Sn::StateLHA; verbose = false)
for i in eachindex(Snplus1.l_var) 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)
...@@ -119,17 +119,12 @@ function read_trajectory(A::LHA, σ::Trajectory; verbose = false) ...@@ -119,17 +119,12 @@ function read_trajectory(A::LHA, σ::Trajectory; verbose = false)
A_new = LHA(A, (σ.m)._map_obs_var_idx) A_new = LHA(A, (σ.m)._map_obs_var_idx)
l_t = times(σ) l_t = times(σ)
l_tr = transitions(σ) l_tr = transitions(σ)
mat_x = zeros(Int, length_states(σ), σ.m.d)
for (i,var) in enumerate(σ.m.g)
mat_x[:,i] = σ[var]
end
Sn = init_state(A_new, σ[1], l_t[1]) Sn = init_state(A_new, σ[1], l_t[1])
Snplus1 = copy(Sn) Snplus1 = copy(Sn)
if verbose println("Init: ") end if verbose println("Init: ") end
if verbose @show Sn end if verbose @show Sn end
for n in 1:length_states(σ) for n in 1:length_states(σ)
xn = view(mat_x, n, :) next_state!(Snplus1, A_new, σ[n], l_t[n], l_tr[n], Sn; verbose = verbose)
next_state!(Snplus1, A_new, xn, l_t[n], l_tr[n], Sn; verbose = verbose)
if Snplus1.loc in A_new.l_loc_final if Snplus1.loc in A_new.l_loc_final
break break
end end
......
core/model.jl
+ 52
− 78
View file @ 146e2afa
...@@ -11,12 +11,18 @@ end ...@@ -11,12 +11,18 @@ end
function _finish_bounded_trajectory!(values::Vector{Vector{Int}}, times::Vector{Float64}, function _finish_bounded_trajectory!(values::Vector{Vector{Int}}, times::Vector{Float64},
transitions::Vector{Transition}, time_bound::Float64) transitions::Vector{Transition}, time_bound::Float64)
for i = eachindex(values) push!(values[i], values[i][end]) end for i = eachindex(values) push!(values[i], values[i][end]) end
push!(times, time_bound) push!(times, time_bound)
push!(transitions, nothing) push!(transitions, nothing)
end end
function _update_values!(values::Vector{Vector{Int}}, times::Vector{Float64}, transitions::Vector{Transition},
xn::Vector{Int}, tn::Float64, tr_n::Transition, idx::Int)
for k = eachindex(values) values[k][idx] = xn[k] end
times[idx] = tn
transitions[idx] = tr_n
end
""" """
`simulate(m)` `simulate(m)`
...@@ -35,7 +41,7 @@ function simulate(m::ContinuousTimeModel) ...@@ -35,7 +41,7 @@ function simulate(m::ContinuousTimeModel)
transitions[1] = nothing transitions[1] = nothing
# Values at time n # Values at time n
n = 1 n = 1
xn = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability xn = m.x0 # View for type stability
tn = m.t0 tn = m.t0
# at time n+1 # at time n+1
isabsorbing::Bool = m.isabsorbing(m.p,xn) isabsorbing::Bool = m.isabsorbing(m.p,xn)
...@@ -48,22 +54,20 @@ function simulate(m::ContinuousTimeModel) ...@@ -48,22 +54,20 @@ function simulate(m::ContinuousTimeModel)
end end
return Trajectory(m, values, times, transitions) return Trajectory(m, values, times, transitions)
end end
# First we fill the allocated array # Alloc of vectors where we stock n+1 values
mat_x = zeros(Int, 1, m.d) vec_x = zeros(Int, m.d)
l_t = Float64[0.0] l_t = Float64[0.0]
l_tr = Transition[nothing] l_tr = Transition[nothing]
# First we fill the allocated array
for i = 2:m.estim_min_states for i = 2:m.estim_min_states
m.f!(mat_x, l_t, l_tr, 1, xn, tn, m.p) m.f!(vec_x, l_t, l_tr, xn, tn, m.p)
tn = l_t[1] tn = l_t[1]
if tn > m.time_bound if tn > m.time_bound
break break
end end
n += 1 n += 1
xn = view(mat_x, 1, :) xn = vec_x
# Updating value _update_values!(full_values, times, transitions, xn, tn, l_tr[1], i)
for k = eachindex(full_values) full_values[k][n] = xn[k] end
times[n] = tn
transitions[n] = l_tr[1]
isabsorbing = m.isabsorbing(m.p,xn) isabsorbing = m.isabsorbing(m.p,xn)
if isabsorbing if isabsorbing
break break
...@@ -79,48 +83,34 @@ function simulate(m::ContinuousTimeModel) ...@@ -79,48 +83,34 @@ function simulate(m::ContinuousTimeModel)
return Trajectory(m, values, times, transitions) return Trajectory(m, values, times, transitions)
end end
# Otherwise, buffering system # Otherwise, buffering system
mat_x = zeros(Int, m.buffer_size, m.d) size_tmp = 0
l_t = zeros(Float64, m.buffer_size)
l_tr = Vector{Transition}(undef, m.buffer_size)
# Alloc buffer values
tmp_full_values = Vector{Vector{Int}}(undef, m.d)
for i = eachindex(tmp_full_values) tmp_full_values[i] = zeros(Int, 0) end
tmp_times = zeros(Float64, 0)
tmp_transitions = Vector{Transition}(undef, 0)
tmp_idx = 0
while !isabsorbing && tn <= m.time_bound while !isabsorbing && tn <= m.time_bound
_resize_trajectory!(tmp_full_values, tmp_times, tmp_transitions, tmp_idx+m.buffer_size) # Alloc buffer
_resize_trajectory!(full_values, times, transitions, m.estim_min_states+size_tmp+m.buffer_size)
i = 0 i = 0
while i < m.buffer_size while i < m.buffer_size
i += 1 i += 1
m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p) m.f!(vec_x, l_t, l_tr, xn, tn, m.p)
tn = l_t[i] tn = l_t[1]
if tn > m.time_bound if tn > m.time_bound
i -= 1 i -= 1
break break
end end
xn = view(mat_x, i, :) xn = vec_x
_update_values!(full_values, times, transitions,
xn, tn, l_tr[1], m.estim_min_states+size_tmp+i)
isabsorbing = m.isabsorbing(m.p,xn) isabsorbing = m.isabsorbing(m.p,xn)
if isabsorbing if isabsorbing
break break
end end
end end
# Update values # If simulation ended before the end of buffer
rng_tmp = (tmp_idx+1):(tmp_idx+m.buffer_size)
for k = eachindex(tmp_full_values) tmp_full_values[k][rng_tmp] = view(mat_x, :, k) end
tmp_times[rng_tmp] = l_t
tmp_transitions[rng_tmp] = l_tr
if i < m.buffer_size if i < m.buffer_size
_resize_trajectory!(tmp_full_values, tmp_times, tmp_transitions, tmp_idx+i) _resize_trajectory!(full_values, times, transitions, m.estim_min_states+size_tmp+i)
end end
tmp_idx += i size_tmp += i
n += i n += i
end end
# Push the temporary values
for k = eachindex(full_values) append!(full_values[k], tmp_full_values[k]) end
append!(times, tmp_times)
append!(transitions, tmp_transitions)
values = full_values[m._g_idx] values = full_values[m._g_idx]
if isbounded(m) if isbounded(m)
# Add last value: the convention is the last transition is nothing, # Add last value: the convention is the last transition is nothing,
...@@ -144,13 +134,13 @@ function simulate(product::SynchronizedModel) ...@@ -144,13 +134,13 @@ function simulate(product::SynchronizedModel)
S0 = init_state(A, m.x0, m.t0) S0 = init_state(A, m.x0, m.t0)
# Values at time n # Values at time n
n = 1 n = 1
xn = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability xn = m.x0 # View for type stability
tn = m.t0 tn = m.t0
Sn = copy(S0) Sn = copy(S0)
# at time n+1
isabsorbing::Bool = m.isabsorbing(m.p,xn) isabsorbing::Bool = m.isabsorbing(m.p,xn)
isacceptedLHA::Bool = isaccepted(Sn) isacceptedLHA::Bool = isaccepted(Sn)
if isabsorbing || isacceptedLHA # If x0 is absorbing
if isabsorbing || isacceptedLHA
_resize_trajectory!(full_values, times, transitions, 1) _resize_trajectory!(full_values, times, transitions, 1)
values = full_values[m._g_idx] values = full_values[m._g_idx]
if isbounded(m) if isbounded(m)
...@@ -158,95 +148,79 @@ function simulate(product::SynchronizedModel) ...@@ -158,95 +148,79 @@ function simulate(product::SynchronizedModel)
end end
return SynchronizedTrajectory(Sn, product, values, times, transitions) return SynchronizedTrajectory(Sn, product, values, times, transitions)
end end
# First we fill the allocated array # Alloc of vectors where we stock n+1 values
mat_x = zeros(Int, 1, m.d) vec_x = zeros(Int, m.d)
l_t = Float64[0.0] l_t = Float64[0.0]
l_tr = Transition[nothing] l_tr = Transition[nothing]
Snplus1 = copy(Sn) Snplus1 = copy(Sn)
# First we fill the allocated array
for i = 2:m.estim_min_states for i = 2:m.estim_min_states
m.f!(mat_x, l_t, l_tr, 1, xn, tn, m.p) m.f!(vec_x, l_t, l_tr, xn, tn, m.p)
tn = l_t[1] tn = l_t[1]
if tn > m.time_bound if tn > m.time_bound
break break
end end
n += 1 n += 1
xn = view(mat_x, 1, :) xn = vec_x
tr_n = l_tr[1] tr_n = l_tr[1]
next_state!(Snplus1, A, xn, tn, tr_n, Sn) next_state!(Snplus1, A, xn, tn, tr_n, Sn)
_update_values!(full_values, times, transitions, xn, tn, tr_n, i)
Sn = Snplus1 Sn = Snplus1
# Updating value
for k = eachindex(full_values) full_values[k][n] = xn[k] end
times[n] = tn
transitions[n] = l_tr[1]
isabsorbing = m.isabsorbing(m.p,xn) isabsorbing = m.isabsorbing(m.p,xn)
isacceptedLHA = isaccepted(Snplus1) isacceptedLHA = isaccepted(Snplus1)
if isabsorbing || isacceptedLHA if isabsorbing || isacceptedLHA
break break
end end
end end
# If simulation ended # If simulation ended before the estimation of states
if n < m.estim_min_states if n < m.estim_min_states
_resize_trajectory!(full_values, times, transitions, n) _resize_trajectory!(full_values, times, transitions, n)
values = full_values[m._g_idx] values = full_values[m._g_idx]
if isbounded(m) if isbounded(m)
_finish_bounded_trajectory!(values, times, transitions, m.time_bound) _finish_bounded_trajectory!(values, times, transitions, m.time_bound)
end end
return SynchronizedTrajectory(Snplus1, product, values, times, transitions) return SynchronizedTrajectory(Sn, product, values, times, transitions)
end end
# Otherwise, buffering system # Otherwise, buffering system
mat_x = zeros(Int, m.buffer_size, m.d) size_tmp = 0
l_t = zeros(Float64, m.buffer_size)
l_tr = Vector{Transition}(undef, m.buffer_size)
# Alloc buffer values
tmp_full_values = Vector{Vector{Int}}(undef, m.d)
for i = eachindex(tmp_full_values) tmp_full_values[i] = zeros(Int, 0) end
tmp_times = zeros(Float64, 0)
tmp_transitions = Vector{Transition}(undef, 0)
tmp_idx = 0
while !isabsorbing && tn <= m.time_bound && !isacceptedLHA while !isabsorbing && tn <= m.time_bound && !isacceptedLHA
_resize_trajectory!(tmp_full_values, tmp_times, tmp_transitions, tmp_idx+m.buffer_size) # Alloc buffer
_resize_trajectory!(full_values, times, transitions, m.estim_min_states+size_tmp+m.buffer_size)
i = 0 i = 0
while i < m.buffer_size while i < m.buffer_size
i += 1 i += 1
m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p) m.f!(vec_x, l_t, l_tr, xn, tn, m.p)
tn = l_t[i] tn = l_t[1]
if tn > m.time_bound if tn > m.time_bound
i -= 1 i -= 1
break break
end end
xn = view(mat_x, i, :) xn = vec_x
tr_n = l_tr[i] tr_n = l_tr[1]
next_state!(Snplus1, A, xn, tn, tr_n, Sn) 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 Sn = Snplus1
isabsorbing = m.isabsorbing(m.p,xn) isabsorbing = m.isabsorbing(m.p,xn)
isacceptedLHA = isaccepted(Snplus1) isacceptedLHA = isaccepted(Snplus1)
if isabsorbing || isacceptedLHA if isabsorbing || isacceptedLHA
break break
end end
end end
# Update values # If simulation ended before the end of buffer
rng_tmp = (tmp_idx+1):(tmp_idx+m.buffer_size)
for k = eachindex(tmp_full_values) tmp_full_values[k][rng_tmp] = view(mat_x, :, k) end
tmp_times[rng_tmp] = l_t
tmp_transitions[rng_tmp] = l_tr
if i < m.buffer_size if i < m.buffer_size
_resize_trajectory!(tmp_full_values, tmp_times, tmp_transitions, tmp_idx+i) _resize_trajectory!(full_values, times, transitions, m.estim_min_states+size_tmp+i)
end end
tmp_idx += i size_tmp += i
n += i n += i
end end
# Push the temporary values
for k = eachindex(full_values) append!(full_values[k], tmp_full_values[k]) end
append!(times, tmp_times)
append!(transitions, tmp_transitions)
values = full_values[m._g_idx] values = full_values[m._g_idx]
if isbounded(m) && !isaccepted(Sn) if isbounded(m) && !isaccepted(Sn)
# Add last value: the convention is the last transition is nothing, # Add last value: the convention is the last transition is nothing,
# the trajectory is bounded # the trajectory is bounded
_finish_bounded_trajectory!(values, times, transitions, m.time_bound) _finish_bounded_trajectory!(values, times, transitions, m.time_bound)
end end
return SynchronizedTrajectory(Snplus1, product, values, times, transitions) return SynchronizedTrajectory(Sn, product, values, times, transitions)
end end
...@@ -274,7 +248,7 @@ function check_consistency(m::ContinuousTimeModel) ...@@ -274,7 +248,7 @@ function check_consistency(m::ContinuousTimeModel)
@assert length(m.g) <= m.d @assert length(m.g) <= m.d
@assert length(m._g_idx) == length(m.g) @assert length(m._g_idx) == length(m.g)
@assert m.buffer_size >= 0 @assert m.buffer_size >= 0
@assert typeof(m.isabsorbing(m.p, view(reshape(m.x0, 1, m.d), 1, :))) == Bool @assert typeof(m.isabsorbing(m.p, m.x0)) == Bool
return true return true
end end
......
models/ER.jl
+ 16
− 13
View file @ 146e2afa
import StaticArrays: SVector, SMatrix, @SMatrix import StaticArrays: SVector, SMatrix, @SVector, @SMatrix
d=4 d=4
k=3 k=3
...@@ -9,18 +9,21 @@ l_tr_ER = ["R1","R2","R3"] ...@@ -9,18 +9,21 @@ l_tr_ER = ["R1","R2","R3"]
p_ER = [1.0, 1.0, 1.0] p_ER = [1.0, 1.0, 1.0]
x0_ER = [100, 100, 0, 0] x0_ER = [100, 100, 0, 0]
t0_ER = 0.0 t0_ER = 0.0
function ER_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}}, idx::Int, function ER_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}},
xn::SubArray{Int,1}, tn::Float64, p::Vector{Float64}) xn::Vector{Int}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2] a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[3] a2 = p[2] * xn[3]
a3 = p[3] * xn[3] a3 = p[3] * xn[3]
l_a = SVector(a1, a2, a3) l_a = SVector(a1, a2, a3)
asum = sum(l_a) asum = sum(l_a)
l_nu = @SMatrix [-1 1 1; nu_1 = SVector(-1, -1, 1, 0)
-1 1 0; nu_2 = SVector(1, 1, -1, 0)
1 -1 -1; nu_3 = SVector(1, 0, -1, 1)
0 0 1] l_nu = SVector(nu_1, nu_2, nu_3)
u1, u2 = rand(), rand() l_str_R = SVector("R1", "R2", "R3")
u1 = rand()
u2 = rand()
tau = - log(u1) / asum tau = - log(u1) / asum
b_inf = 0.0 b_inf = 0.0
b_sup = a1 b_sup = a1
...@@ -34,14 +37,14 @@ function ER_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Noth ...@@ -34,14 +37,14 @@ function ER_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Noth
b_sup += l_a[i+1] b_sup += l_a[i+1]
end end
nu = @view l_nu[:,reaction] # macro for avoiding a copy nu = l_nu[reaction]
for i = 1:4 for i = 1:4
mat_x[idx,i] = xn[i]+nu[i] xnplus1[i] = xn[i]+nu[i]
end end
l_t[idx] = tn + tau l_t[1] = tn + tau
l_tr[idx] = "R$(reaction)" l_tr[1] = l_str_R[reaction]
end end
isabsorbing_ER(p::Vector{Float64},xn::SubArray{Int,1}) = isabsorbing_ER(p::Vector{Float64},xn::Vector{Int}) =
(p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0 (p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0
g_ER = ["P"] g_ER = ["P"]
......
models/SIR.jl
+ 15
− 13
View file @ 146e2afa
import StaticArrays: SVector, SMatrix, @SMatrix import StaticArrays: SVector, SMatrix, @SVector, @SMatrix
d=3 d=3
k=2 k=2
...@@ -9,18 +9,20 @@ l_tr_SIR = ["R1","R2"] ...@@ -9,18 +9,20 @@ l_tr_SIR = ["R1","R2"]
p_SIR = [0.0012, 0.05] p_SIR = [0.0012, 0.05]
x0_SIR = [95, 5, 0] x0_SIR = [95, 5, 0]
t0_SIR = 0.0 t0_SIR = 0.0
function SIR_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}}, idx::Int, function SIR_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}},
xn::SubArray{Int,1}, tn::Float64, p::Vector{Float64}) xn::Vector{Int}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2] a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[2] a2 = p[2] * xn[2]
l_a = SVector(a1, a2) l_a = SVector(a1, a2)
asum = sum(l_a) asum = sum(l_a)
# column-major order # column-major order
l_nu = @SMatrix [-1 0; nu_1 = SVector(-1, 1, 0)
1 -1; nu_2 = SVector(0, -1, 1)
0 1] l_nu = SVector(nu_1, nu_2)
l_str_R = SVector("R1","R2")
u1, u2 = rand(), rand()
u1 = rand()
u2 = rand()
tau = - log(u1) / asum tau = - log(u1) / asum
b_inf = 0.0 b_inf = 0.0
b_sup = a1 b_sup = a1
...@@ -34,14 +36,14 @@ function SIR_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Not ...@@ -34,14 +36,14 @@ function SIR_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Not
b_sup += l_a[i+1] b_sup += l_a[i+1]
end end
nu = view(l_nu, :, reaction) # macro for avoiding a copy nu = l_nu[reaction]
for i = 1:3 for i = 1:3
mat_x[idx,i] = xn[i]+nu[i] xnplus1[i] = xn[i]+nu[i]
end end
l_t[idx] = tn + tau l_t[1] = tn + tau
l_tr[idx] = "R$(reaction)" l_tr[1] = l_str_R[reaction]
end end
isabsorbing_SIR(p::Vector{Float64}, xn::SubArray{Int,1}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0 isabsorbing_SIR(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
g_SIR = ["I"] g_SIR = ["I"]
SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_SIR,p_SIR,x0_SIR,t0_SIR,SIR_f!,isabsorbing_SIR; g=g_SIR) SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_SIR,p_SIR,x0_SIR,t0_SIR,SIR_f!,isabsorbing_SIR; g=g_SIR)
......
tests/cosmos/distance_F/ER_1D.jl
+ 7
− 2
View file @ 146e2afa
...@@ -59,9 +59,14 @@ for exp in l_exp ...@@ -59,9 +59,14 @@ for exp in l_exp
l_dist_i_threads[Threads.threadid()] += (σ.S)["d"] l_dist_i_threads[Threads.threadid()] += (σ.S)["d"]
end end
l_dist_pkg[i] = sum(l_dist_i_threads) / nb_sim l_dist_pkg[i] = sum(l_dist_i_threads) / nb_sim
global test_all = test_all && isapprox(l_dist_cosmos[i], l_dist_pkg[i]; atol = width) 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
end
end end
end end
rm("Result_dist_F_$(str_model).res") rm("Result_dist_F_$(str_model).res")
......
tests/dist_lp/dist_case_2.jl
+ 0
− 1
View file @ 146e2afa
...@@ -35,7 +35,6 @@ for p = 1:2 ...@@ -35,7 +35,6 @@ for p = 1:2
res_int = int^(1/p) res_int = int^(1/p)
test_2 = isapprox(res, res_int; atol = err) test_2 = isapprox(res, res_int; atol = err)
global test_all = test_all && test_1 && test_1_bis && test_2 global test_all = test_all && test_1 && test_1_bis && test_2
end end
......
tests/dist_lp/dist_l1_case_1.jl
+ 38
− 32
View file @ 146e2afa
...@@ -5,36 +5,42 @@ load_model("SIR") ...@@ -5,36 +5,42 @@ load_model("SIR")
# Case 1 : 6.4 # Case 1 : 6.4
x_obs = [2, 4, 3, 3] test_all = true
t_x = [0.0, 0.5, 1.2, 2.2]
t_y_bis = [0.0, 0.5, 1.2, 2.2] let x_obs, y_obs, t_x, t_y, σ1, σ2
values = [zeros(length(x_obs))] x_obs = [2, 4, 3, 3]
values[1] = x_obs t_x = [0.0, 0.5, 1.2, 2.2]
l_tr = Vector{Nothing}(nothing, length(x_obs)) t_y_bis = [0.0, 0.5, 1.2, 2.2]
σ1 = Trajectory(SIR, values, t_x, l_tr) values = [zeros(length(x_obs))]
values[1] = x_obs
y_obs = [6, 6] l_tr = Vector{Nothing}(nothing, length(x_obs))
t_y = [0.0, 2.2] σ1 = Trajectory(SIR, values, t_x, l_tr)
t_x_bis = [0.0, 2.2]
values = [zeros(length(y_obs))] y_obs = [6, 6]
values[1] = y_obs t_y = [0.0, 2.2]
l_tr = Vector{Nothing}(nothing, length(y_obs)) t_x_bis = [0.0, 2.2]
σ2 = Trajectory(SIR, values, t_y, l_tr) values = [zeros(length(y_obs))]
values[1] = y_obs
test_1 = dist_lp(x_obs, t_x, y_obs, t_y; p=1) == 6.4 l_tr = Vector{Nothing}(nothing, length(y_obs))
test_1 = test_1 && dist_lp(σ1, σ2; p=1) == 6.4 σ2 = Trajectory(SIR, values, t_y, l_tr)
f_x(t::Real) = MarkovProcesses._f_step(x_obs, t_x, t) test_1 = dist_lp(x_obs, t_x, y_obs, t_y; p=1) == 6.4
f_y(t::Real) = MarkovProcesses._f_step(y_obs, t_y, t) test_1 = test_1 && dist_lp(σ1, σ2; p=1) == 6.4
diff_f(t) = abs(f_x(t) - f_y(t))
int, err = quadgk(diff_f, 0.0, 2.2) f_x(t::Real) = MarkovProcesses._f_step(x_obs, t_x, t)
f_y(t::Real) = MarkovProcesses._f_step(y_obs, t_y, t)
test_2 = isapprox(6.4, int; atol=err) diff_f(t) = abs(f_x(t) - f_y(t))
int, err = quadgk(diff_f, 0.0, 2.2)
# Case 1 bis : inverse of case 1
test_2 = isapprox(6.4, int; atol=err)
test_1_bis = dist_lp(y_obs, t_y, x_obs, t_x; p=1) == 6.4
test_1_bis = test_1_bis && dist_lp(σ2, σ1; p=1) == 6.4 # Case 1 bis : inverse of case 1
return test_1 && test_1_bis && test_2 test_1_bis = dist_lp(y_obs, t_y, x_obs, t_x; p=1) == 6.4
test_1_bis = test_1_bis && dist_lp(σ2, σ1; p=1) == 6.4
global test_all = test_1 && test_1_bis && test_2
end
return test_all
tests/unit/check_trajectory_consistency.jl
+ 6
− 6
View file @ 146e2afa
...@@ -24,8 +24,8 @@ end ...@@ -24,8 +24,8 @@ end
σ2_dump = simulate(ER) σ2_dump = simulate(ER)
for i = 1:nb_sim for i = 1:nb_sim
σ = simulate(SIR) local σ = simulate(SIR)
σ2 = simulate(ER) local σ2 = simulate(ER)
try try
global test_all = test_all && check_consistency(σ) && check_consistency(σ2) && global test_all = test_all && check_consistency(σ) && check_consistency(σ2) &&
!isbounded(σ) && !isbounded(σ2) !isbounded(σ) && !isbounded(σ2)
...@@ -48,8 +48,8 @@ end ...@@ -48,8 +48,8 @@ end
SIR.time_bound = 1.0 SIR.time_bound = 1.0
ER.time_bound = 0.01 ER.time_bound = 0.01
for i = 1:nb_sim for i = 1:nb_sim
σ = simulate(SIR) local σ = simulate(SIR)
σ2 = simulate(ER) local σ2 = simulate(ER)
try try
global test_all = test_all && check_consistency(σ) && check_consistency(σ2) && global test_all = test_all && check_consistency(σ) && check_consistency(σ2) &&
isbounded(σ) && isbounded(σ2) isbounded(σ) && isbounded(σ2)
...@@ -71,8 +71,8 @@ end ...@@ -71,8 +71,8 @@ end
SIR.time_bound = Inf SIR.time_bound = Inf
ER.time_bound = Inf ER.time_bound = Inf
new_x0_SIR = view(reshape([95, 0, 5], 1, :), 1, :) new_x0_SIR = [95, 0, 5]
new_x0_ER = view(reshape([0, 0, 0, 100], 1, :), 1, :) new_x0_ER = [0, 0, 0, 100]
SIR.x0 = new_x0_SIR SIR.x0 = new_x0_SIR
ER.x0 = new_x0_ER ER.x0 = new_x0_ER
σ = simulate(SIR) σ = simulate(SIR)
......
tests/unit/is_always_bounded_sir.jl
+ 1
− 1
View file @ 146e2afa
...@@ -5,7 +5,7 @@ test = true ...@@ -5,7 +5,7 @@ test = true
load_model("SIR") load_model("SIR")
SIR.time_bound = 120.0 SIR.time_bound = 120.0
for i = 1:1000 for i = 1:1000
σ = simulate(SIR) local σ = simulate(SIR)
global test = test && isbounded(σ) global test = test && isbounded(σ)
end end
......
tests/unit/long_sim_er.jl 0 → 100644
+ 8
− 0
View file @ 146e2afa
using MarkovProcesses
load_model("ER")
set_param!(ER, "k1", 0.2)
set_param!(ER, "k2", 40.0)
return true
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment