From 40cf91a44e48fd14c2198c8c73cfb980a68f4ccb Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Sun, 7 Feb 2021 15:45:05 +0100
Subject: [PATCH] Several fixs: integrated poisson model wasn't updated
improvement of plots fix of synchronized simulation when state is absorbing
All tests passed
---
core/model.jl | 22 ++++++++++++++--------
core/plots.jl | 12 +++++++-----
models/poisson.jl | 8 ++++----
tests/automata/absorbing_x0_p.jl | 11 +++++++++--
4 files changed, 34 insertions(+), 19 deletions(-)
diff --git a/core/model.jl b/core/model.jl
index 4e2a797..83b5313 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -133,7 +133,10 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
if p != nothing
p_sim = p
end
+ x0 = getfield(m, :x0)
+ t0 = getfield(m, :t0)
time_bound = getfield(m, :time_bound)
+ S0 = init_state(A, x0, t0)
buffer_size = getfield(m, :buffer_size)
estim_min_states = getfield(m, :estim_min_states)
# First alloc
@@ -142,15 +145,15 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl
times = zeros(Float64, estim_min_states)
transitions = Vector{Transition}(undef, estim_min_states)
# Initial values
- for i = eachindex(full_values) full_values[i][1] = getfield(m, :x0)[i] end
- times[1] = getfield(m, :t0)
+ for i = eachindex(full_values) full_values[i][1] = x0[i] end
+ times[1] = t0
transitions[1] = nothing
- S0 = init_state(A, getfield(m, :x0), getfield(m, :t0))
# Values at time n
n = 1
- xn = copy(getfield(m, :x0))
- tn = getfield(m, :t0)
+ xn = copy(x0)
+ tn = copy(t0)
Sn = copy(S0)
+ next_state!(Sn, A, x0, t0, nothing, S0, x0, p_sim; verbose = verbose)
isabsorbing::Bool = getfield(m, :isabsorbing)(p_sim,xn)
isacceptedLHA::Bool = isaccepted(Sn)
# Alloc of vectors where we stock n+1 values
@@ -269,13 +272,16 @@ function volatile_simulate(product::SynchronizedModel;
if p != nothing
p_sim = p
end
+ x0 = getfield(m, :x0)
+ t0 = getfield(m, :t0)
time_bound = getfield(m, :time_bound)
- S0 = init_state(A, getfield(m, :x0), getfield(m, :t0))
+ S0 = init_state(A, x0, t0)
# Values at time n
n = 1
- xn = copy(getfield(m, :x0))
- tn = getfield(m, :t0)
+ xn = copy(x0)
+ tn = copy(t0)
Sn = copy(S0)
+ next_state!(Sn, A, x0, t0, nothing, S0, x0, p_sim; verbose = verbose)
isabsorbing::Bool = getfield(m, :isabsorbing)(p_sim,xn)
isacceptedLHA::Bool = isaccepted(Sn)
# Alloc of vectors where we stock n+1 values
diff --git a/core/plots.jl b/core/plots.jl
index 2f260ff..7e1a15f 100644
--- a/core/plots.jl
+++ b/core/plots.jl
@@ -53,25 +53,27 @@ function plot(σ::AbstractTrajectory, vars::VariableModel...; plot_transitions::
end
end
-function plot!(A::LHA)
+function plot!(A::LHA; label::String = "")
+ label_region = (label == "") ? "Region LHA" : label
if A.name in ["F property", "G property"]
if haskey(A.constants, :x1) && haskey(A.constants, :x2) && haskey(A.constants, :t1) && haskey(A.constants, :t2)
x1, x2, t1, t2 = A.constants[:x1], A.constants[:x2], A.constants[:t1], A.constants[:t2]
- plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5, label = "Region LHA")
+ plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5, label = label_region)
end
end
if A.name in ["G and F property"]
if haskey(A.constants, :x1) && haskey(A.constants, :x2) && haskey(A.constants, :t1) && haskey(A.constants, :t2)
x1, x2, t1, t2 = A.constants[:x1], A.constants[:x2], A.constants[:t1], A.constants[:t2]
- plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5, label = "Region LHA")
+ plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5, label = label_region)
end
if haskey(A.constants, :x3) && haskey(A.constants, :x4) && haskey(A.constants, :t3) && haskey(A.constants, :t4)
x3, x4, t3, t4 = A.constants[:x3], A.constants[:x4], A.constants[:t3], A.constants[:t4]
- plot!(Shape([(t3,x3), (t3,x4), (t4,x4), (t4,x3), (t3,x3)]), opacity = 0.5, label = "Region LHA")
+ plot!(Shape([(t3,x3), (t3,x4), (t4,x4), (t4,x3), (t3,x3)]), opacity = 0.5, label = label_region)
end
end
+ label_region = (label == "") ? "L, H" : label
if A.name == "Period"
- hline!([A.constants[:L], A.constants[:H]], label = "L, H", linestyle = :dot)
+ hline!([A.constants[:L], A.constants[:H]], label = label_region, linestyle = :dot)
end
end
diff --git a/models/poisson.jl b/models/poisson.jl
index c0133ae..02de095 100644
--- a/models/poisson.jl
+++ b/models/poisson.jl
@@ -10,8 +10,8 @@ l_tr_poisson = [:R]
p_poisson = [5.0]
x0_poisson = [0]
t0_poisson = 0.0
-function poisson_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition},
- xn::Vector{Int}, tn::Float64, p::Vector{Float64})
+@everywhere function poisson_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition},
+ xn::Vector{Int}, tn::Float64, p::Vector{Float64})
u1 = rand()
tau = (-log(u1)/p[1])
@@ -19,12 +19,12 @@ function poisson_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Tra
l_t[1] = tn + tau
l_tr[1] = :R
end
-isabsorbing_poisson(p::Vector{Float64}, xn::Vector{Int}) = p[1] === 0.0
+@everywhere isabsorbing_poisson(p::Vector{Float64}, xn::Vector{Int}) = p[1] === 0.0
g_poisson = [:N]
poisson = ContinuousTimeModel(d,k,dict_var_poisson,dict_p_poisson,l_tr_poisson,
p_poisson,x0_poisson,t0_poisson,
- poisson_f!,isabsorbing_poisson;
+ getfield(Main, :poisson_f!), getfield(Main, :isabsorbing_poisson);
g=g_poisson, time_bound=1.0, name="Poisson process pkg")
function create_poisson(new_p::Vector{Float64})
poisson_new = deepcopy(poisson)
diff --git a/tests/automata/absorbing_x0_p.jl b/tests/automata/absorbing_x0_p.jl
index 2fd9843..14c6219 100644
--- a/tests/automata/absorbing_x0_p.jl
+++ b/tests/automata/absorbing_x0_p.jl
@@ -1,6 +1,8 @@
using MarkovProcesses
+test_all = true
+
load_model("ER")
observe_all!(ER)
load_automaton("automaton_F")
@@ -14,7 +16,6 @@ A_G_F_R6 = create_automaton_G_and_F(ER, x1, x2, t1, t2, :E,
x3, x4, t3, t4, :P)
set_param!(ER, [0.0,0.0,0.0])
-test_all = true
for A in [A_F_R1, A_G_R5, A_G_F_R6]
local sync = A*ER
local σ, σ2 = simulate(sync), simulate(sync)
@@ -23,12 +24,18 @@ end
set_param!(ER, [1.0,1.0,1.0])
set_x0!(ER, [0, 0, 0, 5])
-test_all = true
for A in [A_F_R1, A_G_R5, A_G_F_R6]
local sync = A*ER
local σ, σ2 = simulate(sync), simulate(sync)
global test_all = test_all && isaccepted(σ) && isaccepted(σ2)
end
+load_model("poisson")
+set_param!(poisson, [0.0])
+aut = create_automaton_F(poisson, 5.0, 7.0, 0.75, 1.0, :N)
+sync_poisson = aut * poisson
+test_all = test_all && simulate(sync_poisson).state_lha_end[:d] > 0
+test_all = test_all && volatile_simulate(sync_poisson)[:d] > 0
+
return true
--
GitLab