diff --git a/automata/automaton_F.jl b/automata/automaton_F.jl
index c72c6bc112b16cc21e08071b29dd36486d96bac0..d5d452a3a1a7bace4942bfaee0f734b9eeea9d79 100644
--- a/automata/automaton_F.jl
+++ b/automata/automaton_F.jl
@@ -32,7 +32,7 @@ function create_automaton_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
# "cc" as check_constraints
tuple = ("l0", "l1")
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))
edge1 = Edge([nothing], cc_aut_F_l0l1_1, us_aut_F_l0l1_1!)
map_edges[tuple] = [edge1]
@@ -41,34 +41,34 @@ function create_automaton_F(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
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::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!)
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!)
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!)
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::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!)
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"])
- 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),
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!)
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"])
- 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")
edge3 = Edge([nothing], cc_aut_F_l1l3_3, us_aut_F_l1l3_3!)
map_edges[tuple] = [edge1, edge2, edge3]
@@ -76,11 +76,11 @@ 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::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!)
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::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!)
map_edges[tuple] = [edge1, edge2]
diff --git a/automata/automaton_G.jl b/automata/automaton_G.jl
index 1607b598135df6b95f32f2ae3cf9181de80ab1ca..38f9ded0307b4a3789e35fe3abddf891463f051b 100644
--- a/automata/automaton_G.jl
+++ b/automata/automaton_G.jl
@@ -33,7 +33,7 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
# 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::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!)
map_edges[tuple] = [edge1]
@@ -41,64 +41,64 @@ function create_automaton_G(m::ContinuousTimeModel, x1::Float64, x2::Float64, t1
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::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)
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::SubArray{Int,1}) =
+ 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::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!)
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::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!)
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::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!)
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::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!)
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::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!)
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!)
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::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!)
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::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!)
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!)
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::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["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!)
@@ -106,7 +106,7 @@ 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"])
- 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)
edge1 = Edge([nothing], cc_aut_G_l4l2_1, us_aut_G_l4l2_1!)
map_edges[tuple] = [edge1]
diff --git a/bench/pygmalion/multiple_long_sim_er.jl b/bench/pygmalion/multiple_long_sim_er.jl
index f804a69355db94b91d6efb297b882587fdcff7fd..98254998cccd5abc37d1e318db810748043390ad 100644
--- a/bench/pygmalion/multiple_long_sim_er.jl
+++ b/bench/pygmalion/multiple_long_sim_er.jl
@@ -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)
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(b2_pyg), mean(b2_pyg), maximum(b2_pyg)
+#@show minimum(b2_pyg), mean(b2_pyg), maximum(b2_pyg)
println("MarkovProcesses:")
using MarkovProcesses
MarkovProcesses.load_model("ER")
ER.time_bound = 20.0
-ER.estim_min_states = 9000
set_param!(ER, "k1", 0.2)
set_param!(ER, "k2", 40.0)
@timev MarkovProcesses.simulate(ER)
println("Default buffer size=10")
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(b2), mean(b2), maximum(b2)
-println("Buffer size 100")
+#@show minimum(b2), mean(b2), maximum(b2)
+println("Buffer size 100 / min_states 8000")
+ER.estim_min_states = 8000
ER.buffer_size = 100
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(b2_100), mean(b2_100), maximum(b2_100)
+#@show minimum(b2_100), mean(b2_100), maximum(b2_100)
diff --git a/core/_tests_simulate.jl b/core/_tests_simulate.jl
index cdc1ac7d6a58b55284b2b54a4d453756caf27e8c..b5bd48051690186bb88827d3098298cbf2459816 100644
--- a/core/_tests_simulate.jl
+++ b/core/_tests_simulate.jl
@@ -334,66 +334,3 @@ function _simulate_d7458(m::ContinuousTimeModel)
return Trajectory(m, values, times, transitions)
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
-
-
diff --git a/core/lha.jl b/core/lha.jl
index 2445b23d076bbcb53096adf10f213ba725ee0cd0..9e146da79f779f348ee1270fa14846729ef7d23e 100644
--- a/core/lha.jl
+++ b/core/lha.jl
@@ -2,7 +2,7 @@
load_automaton(automaton::String) = include(get_module_path() * "/automata/$(automaton).jl")
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)
# 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)
end
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)
for i in eachindex(Snplus1.l_var)
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)
A_new = LHA(A, (σ.m)._map_obs_var_idx)
l_t = times(σ)
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])
Snplus1 = copy(Sn)
if verbose println("Init: ") end
if verbose @show Sn end
for n in 1:length_states(σ)
- xn = view(mat_x, n, :)
- next_state!(Snplus1, A_new, xn, l_t[n], l_tr[n], Sn; verbose = verbose)
+ next_state!(Snplus1, A_new, σ[n], l_t[n], l_tr[n], Sn; verbose = verbose)
if Snplus1.loc in A_new.l_loc_final
break
end
diff --git a/core/model.jl b/core/model.jl
index 7865a8b44e4b4cb2563131a8378285e196309cf4..3e363ca8603f730ddfe91ce46ed7181166f0ce1e 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -11,12 +11,18 @@ end
function _finish_bounded_trajectory!(values::Vector{Vector{Int}}, times::Vector{Float64},
transitions::Vector{Transition}, time_bound::Float64)
-
for i = eachindex(values) push!(values[i], values[i][end]) end
push!(times, time_bound)
push!(transitions, nothing)
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)`
@@ -35,7 +41,7 @@ function simulate(m::ContinuousTimeModel)
transitions[1] = nothing
# Values at time n
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
# at time n+1
isabsorbing::Bool = m.isabsorbing(m.p,xn)
@@ -48,22 +54,20 @@ function simulate(m::ContinuousTimeModel)
end
return Trajectory(m, values, times, transitions)
end
- # First we fill the allocated array
- mat_x = zeros(Int, 1, m.d)
+ # Alloc of vectors where we stock n+1 values
+ vec_x = zeros(Int, m.d)
l_t = Float64[0.0]
l_tr = Transition[nothing]
+ # First we fill the allocated array
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]
if tn > m.time_bound
break
end
n += 1
- xn = view(mat_x, 1, :)
- # Updating value
- for k = eachindex(full_values) full_values[k][n] = xn[k] end
- times[n] = tn
- transitions[n] = l_tr[1]
+ xn = vec_x
+ _update_values!(full_values, times, transitions, xn, tn, l_tr[1], i)
isabsorbing = m.isabsorbing(m.p,xn)
if isabsorbing
break
@@ -79,48 +83,34 @@ function simulate(m::ContinuousTimeModel)
return Trajectory(m, values, times, transitions)
end
# Otherwise, buffering system
- mat_x = zeros(Int, m.buffer_size, m.d)
- 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
+ size_tmp = 0
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
while i < m.buffer_size
i += 1
- m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
- tn = l_t[i]
+ m.f!(vec_x, l_t, l_tr, xn, tn, m.p)
+ tn = l_t[1]
if tn > m.time_bound
i -= 1
break
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)
if isabsorbing
break
end
end
- # Update values
- 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 simulation ended before the end of buffer
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
- tmp_idx += i
+ size_tmp += i
n += i
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]
if isbounded(m)
# Add last value: the convention is the last transition is nothing,
@@ -144,13 +134,13 @@ function simulate(product::SynchronizedModel)
S0 = init_state(A, m.x0, m.t0)
# Values at time n
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
Sn = copy(S0)
- # at time n+1
isabsorbing::Bool = m.isabsorbing(m.p,xn)
isacceptedLHA::Bool = isaccepted(Sn)
- if isabsorbing || isacceptedLHA
+ # If x0 is absorbing
+ if isabsorbing || isacceptedLHA
_resize_trajectory!(full_values, times, transitions, 1)
values = full_values[m._g_idx]
if isbounded(m)
@@ -158,95 +148,79 @@ function simulate(product::SynchronizedModel)
end
return SynchronizedTrajectory(Sn, product, values, times, transitions)
end
- # First we fill the allocated array
- mat_x = zeros(Int, 1, m.d)
+ # 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!(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]
if tn > m.time_bound
break
end
n += 1
- xn = view(mat_x, 1, :)
+ 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
- # 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)
isacceptedLHA = isaccepted(Snplus1)
if isabsorbing ||Â isacceptedLHA
break
end
end
- # If simulation ended
+ # If simulation ended before the estimation of states
if n < m.estim_min_states
_resize_trajectory!(full_values, times, transitions, n)
values = full_values[m._g_idx]
if isbounded(m)
_finish_bounded_trajectory!(values, times, transitions, m.time_bound)
end
- return SynchronizedTrajectory(Snplus1, product, values, times, transitions)
+ return SynchronizedTrajectory(Sn, product, values, times, transitions)
end
# Otherwise, buffering system
- mat_x = zeros(Int, m.buffer_size, m.d)
- 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
+ size_tmp = 0
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
while i < m.buffer_size
i += 1
- m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
- tn = l_t[i]
+ m.f!(vec_x, l_t, l_tr, xn, tn, m.p)
+ tn = l_t[1]
if tn > m.time_bound
i -= 1
break
end
- xn = view(mat_x, i, :)
- tr_n = l_tr[i]
+ 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
isabsorbing = m.isabsorbing(m.p,xn)
isacceptedLHA = isaccepted(Snplus1)
- if isabsorbing ||Â isacceptedLHA
+ if isabsorbing || isacceptedLHA
break
end
end
- # Update values
- 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 simulation ended before the end of buffer
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
- tmp_idx += i
+ size_tmp += i
n += i
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]
if isbounded(m) && !isaccepted(Sn)
# Add last value: the convention is the last transition is nothing,
# the trajectory is bounded
_finish_bounded_trajectory!(values, times, transitions, m.time_bound)
end
- return SynchronizedTrajectory(Snplus1, product, values, times, transitions)
+ return SynchronizedTrajectory(Sn, product, values, times, transitions)
end
@@ -274,7 +248,7 @@ function check_consistency(m::ContinuousTimeModel)
@assert length(m.g) <= m.d
@assert length(m._g_idx) == length(m.g)
@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
end
diff --git a/models/ER.jl b/models/ER.jl
index a376d6de13685f5555380cee2e08f3e513530c49..da45ebaecaf5af5d79d0e7f0f232ad7fffd5a971 100644
--- a/models/ER.jl
+++ b/models/ER.jl
@@ -1,5 +1,5 @@
-import StaticArrays: SVector, SMatrix, @SMatrix
+import StaticArrays: SVector, SMatrix, @SVector, @SMatrix
d=4
k=3
@@ -9,18 +9,21 @@ l_tr_ER = ["R1","R2","R3"]
p_ER = [1.0, 1.0, 1.0]
x0_ER = [100, 100, 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,
- xn::SubArray{Int,1}, tn::Float64, p::Vector{Float64})
+function ER_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}},
+ xn::Vector{Int}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[3]
a3 = p[3] * xn[3]
l_a = SVector(a1, a2, a3)
asum = sum(l_a)
- l_nu = @SMatrix [-1 1 1;
- -1 1 0;
- 1 -1 -1;
- 0 0 1]
- u1, u2 = rand(), rand()
+ nu_1 = SVector(-1, -1, 1, 0)
+ nu_2 = SVector(1, 1, -1, 0)
+ nu_3 = SVector(1, 0, -1, 1)
+ l_nu = SVector(nu_1, nu_2, nu_3)
+ l_str_R = SVector("R1", "R2", "R3")
+
+ u1 = rand()
+ u2 = rand()
tau = - log(u1) / asum
b_inf = 0.0
b_sup = a1
@@ -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]
end
- nu = @view l_nu[:,reaction] # macro for avoiding a copy
+ nu = l_nu[reaction]
for i = 1:4
- mat_x[idx,i] = xn[i]+nu[i]
+ xnplus1[i] = xn[i]+nu[i]
end
- l_t[idx] = tn + tau
- l_tr[idx] = "R$(reaction)"
+ l_t[1] = tn + tau
+ l_tr[1] = l_str_R[reaction]
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
g_ER = ["P"]
diff --git a/models/SIR.jl b/models/SIR.jl
index 48e9354db534cf52c94f26a2f68662c42d171c46..5da6d2c1e53f399e5e32b0d654aea7a30b479e5b 100644
--- a/models/SIR.jl
+++ b/models/SIR.jl
@@ -1,5 +1,5 @@
-import StaticArrays: SVector, SMatrix, @SMatrix
+import StaticArrays: SVector, SMatrix, @SVector, @SMatrix
d=3
k=2
@@ -9,18 +9,20 @@ l_tr_SIR = ["R1","R2"]
p_SIR = [0.0012, 0.05]
x0_SIR = [95, 5, 0]
t0_SIR = 0.0
-function SIR_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}}, idx::Int,
- xn::SubArray{Int,1}, tn::Float64, p::Vector{Float64})
+function SIR_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}},
+ xn::Vector{Int}, tn::Float64, p::Vector{Float64})
a1 = p[1] * xn[1] * xn[2]
a2 = p[2] * xn[2]
l_a = SVector(a1, a2)
asum = sum(l_a)
# column-major order
- l_nu = @SMatrix [-1 0;
- 1 -1;
- 0 1]
-
- u1, u2 = rand(), rand()
+ nu_1 = SVector(-1, 1, 0)
+ nu_2 = SVector(0, -1, 1)
+ l_nu = SVector(nu_1, nu_2)
+ l_str_R = SVector("R1","R2")
+
+ u1 = rand()
+ u2 = rand()
tau = - log(u1) / asum
b_inf = 0.0
b_sup = a1
@@ -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]
end
- nu = view(l_nu, :, reaction) # macro for avoiding a copy
+ nu = l_nu[reaction]
for i = 1:3
- mat_x[idx,i] = xn[i]+nu[i]
+ xnplus1[i] = xn[i]+nu[i]
end
- l_t[idx] = tn + tau
- l_tr[idx] = "R$(reaction)"
+ l_t[1] = tn + tau
+ l_tr[1] = l_str_R[reaction]
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"]
SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_SIR,p_SIR,x0_SIR,t0_SIR,SIR_f!,isabsorbing_SIR; g=g_SIR)
diff --git a/tests/cosmos/distance_F/ER_1D.jl b/tests/cosmos/distance_F/ER_1D.jl
index 96ce2bf6e61226b050c27931a51ab2ebfab95465..0b2cc8c90d5d597f364c68b9b20cf69acec9dabc 100644
--- a/tests/cosmos/distance_F/ER_1D.jl
+++ b/tests/cosmos/distance_F/ER_1D.jl
@@ -59,9 +59,14 @@ for exp in l_exp
l_dist_i_threads[Threads.threadid()] += (σ.S)["d"]
end
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
rm("Result_dist_F_$(str_model).res")
diff --git a/tests/dist_lp/dist_case_2.jl b/tests/dist_lp/dist_case_2.jl
index c66ce3b1dfc4dd66b67339fb312cc24f7eac78be..102280bdcec0602ba10252295a0ceec16df576b8 100644
--- a/tests/dist_lp/dist_case_2.jl
+++ b/tests/dist_lp/dist_case_2.jl
@@ -35,7 +35,6 @@ for p = 1:2
res_int = int^(1/p)
test_2 = isapprox(res, res_int; atol = err)
-
global test_all = test_all && test_1 && test_1_bis && test_2
end
diff --git a/tests/dist_lp/dist_l1_case_1.jl b/tests/dist_lp/dist_l1_case_1.jl
index 262c8c77c0140bbf8e23d71540f73fdbf4a412ab..a4fd4e50a572e49f44564fcd0701c7954d50fefe 100644
--- a/tests/dist_lp/dist_l1_case_1.jl
+++ b/tests/dist_lp/dist_l1_case_1.jl
@@ -5,36 +5,42 @@ load_model("SIR")
# Case 1 : 6.4
-x_obs = [2, 4, 3, 3]
-t_x = [0.0, 0.5, 1.2, 2.2]
-t_y_bis = [0.0, 0.5, 1.2, 2.2]
-values = [zeros(length(x_obs))]
-values[1] = x_obs
-l_tr = Vector{Nothing}(nothing, length(x_obs))
-σ1 = Trajectory(SIR, values, t_x, l_tr)
-
-y_obs = [6, 6]
-t_y = [0.0, 2.2]
-t_x_bis = [0.0, 2.2]
-values = [zeros(length(y_obs))]
-values[1] = y_obs
-l_tr = Vector{Nothing}(nothing, length(y_obs))
-σ2 = Trajectory(SIR, values, t_y, l_tr)
-
-test_1 = dist_lp(x_obs, t_x, y_obs, t_y; p=1) == 6.4
-test_1 = test_1 && dist_lp(σ1, σ2; p=1) == 6.4
-
-f_x(t::Real) = MarkovProcesses._f_step(x_obs, t_x, t)
-f_y(t::Real) = MarkovProcesses._f_step(y_obs, t_y, t)
-diff_f(t) = abs(f_x(t) - f_y(t))
-int, err = quadgk(diff_f, 0.0, 2.2)
-
-test_2 = isapprox(6.4, int; atol=err)
-
-# Case 1 bis : inverse of case 1
-
-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
-
-return test_1 && test_1_bis && test_2
+test_all = true
+
+let x_obs, y_obs, t_x, t_y, σ1, σ2
+ x_obs = [2, 4, 3, 3]
+ t_x = [0.0, 0.5, 1.2, 2.2]
+ t_y_bis = [0.0, 0.5, 1.2, 2.2]
+ values = [zeros(length(x_obs))]
+ values[1] = x_obs
+ l_tr = Vector{Nothing}(nothing, length(x_obs))
+ σ1 = Trajectory(SIR, values, t_x, l_tr)
+
+ y_obs = [6, 6]
+ t_y = [0.0, 2.2]
+ t_x_bis = [0.0, 2.2]
+ values = [zeros(length(y_obs))]
+ values[1] = y_obs
+ l_tr = Vector{Nothing}(nothing, length(y_obs))
+ σ2 = Trajectory(SIR, values, t_y, l_tr)
+
+ test_1 = dist_lp(x_obs, t_x, y_obs, t_y; p=1) == 6.4
+ test_1 = test_1 && dist_lp(σ1, σ2; p=1) == 6.4
+
+ f_x(t::Real) = MarkovProcesses._f_step(x_obs, t_x, t)
+ f_y(t::Real) = MarkovProcesses._f_step(y_obs, t_y, t)
+ diff_f(t) = abs(f_x(t) - f_y(t))
+ int, err = quadgk(diff_f, 0.0, 2.2)
+
+ test_2 = isapprox(6.4, int; atol=err)
+
+ # Case 1 bis : inverse of case 1
+
+ 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
diff --git a/tests/unit/check_trajectory_consistency.jl b/tests/unit/check_trajectory_consistency.jl
index f60a07a0ae4b77a031ffe750aaa7d59e8573ce8c..1b3f39a7a9d5a08f29d18f9e1fdf8a4f1fe51bce 100644
--- a/tests/unit/check_trajectory_consistency.jl
+++ b/tests/unit/check_trajectory_consistency.jl
@@ -24,8 +24,8 @@ end
σ2_dump = simulate(ER)
for i = 1:nb_sim
- σ = simulate(SIR)
- σ2 = simulate(ER)
+ local σ = simulate(SIR)
+ local σ2 = simulate(ER)
try
global test_all = test_all && check_consistency(σ) && check_consistency(σ2) &&
!isbounded(σ) && !isbounded(σ2)
@@ -48,8 +48,8 @@ end
SIR.time_bound = 1.0
ER.time_bound = 0.01
for i = 1:nb_sim
- σ = simulate(SIR)
- σ2 = simulate(ER)
+ local σ = simulate(SIR)
+ local σ2 = simulate(ER)
try
global test_all = test_all && check_consistency(σ) && check_consistency(σ2) &&
isbounded(σ) && isbounded(σ2)
@@ -71,8 +71,8 @@ end
SIR.time_bound = Inf
ER.time_bound = Inf
-new_x0_SIR = view(reshape([95, 0, 5], 1, :), 1, :)
-new_x0_ER = view(reshape([0, 0, 0, 100], 1, :), 1, :)
+new_x0_SIR = [95, 0, 5]
+new_x0_ER = [0, 0, 0, 100]
SIR.x0 = new_x0_SIR
ER.x0 = new_x0_ER
σ = simulate(SIR)
diff --git a/tests/unit/is_always_bounded_sir.jl b/tests/unit/is_always_bounded_sir.jl
index 7b478e061713b922bd44855e05461a7cb00f8e8c..778b8fc0673fbb76a414f85a50e03083500a8b8b 100644
--- a/tests/unit/is_always_bounded_sir.jl
+++ b/tests/unit/is_always_bounded_sir.jl
@@ -5,7 +5,7 @@ test = true
load_model("SIR")
SIR.time_bound = 120.0
for i = 1:1000
- σ = simulate(SIR)
+ local σ = simulate(SIR)
global test = test && isbounded(σ)
end
diff --git a/tests/unit/long_sim_er.jl b/tests/unit/long_sim_er.jl
new file mode 100644
index 0000000000000000000000000000000000000000..ea5342ddd6efbd0578a1b055d6b3503bfe356295
--- /dev/null
+++ b/tests/unit/long_sim_er.jl
@@ -0,0 +1,8 @@
+
+using MarkovProcesses
+load_model("ER")
+set_param!(ER, "k1", 0.2)
+set_param!(ER, "k2", 40.0)
+
+return true
+