diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index 9a17dd2f80af7538a4cd538fc943830815167080..975cb6727200eae29eff6b0f29bb8e3ee89421d5 100644
--- a/core/MarkovProcesses.jl
+++ b/core/MarkovProcesses.jl
@@ -38,7 +38,7 @@ export init_state, next_state!, read_trajectory
export get_index, get_value, length_var, isaccepted
# Model related methods
-export simulate, volatile_simulate
+export simulate, volatile_simulate, change_simulation_stop_criteria
export distribute_mean_value_lha, mean_value_lha, distribute_prob_accept_lha, probability_var_value_lha
export number_simulations_smc_chernoff, smc_chernoff
export set_param!, set_x0!, set_time_bound!, set_observed_var!, observe_all!
diff --git a/core/model.jl b/core/model.jl
index 9f6199f4b483075afd90818ddb963afe1f2f21ce..7b8c9153810f8539bfc95930142284ce166bc7f9 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -271,12 +271,12 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym
n += i
end
values = full_values[getfield(m, :_g_idx)]
- if isbounded(m) && !isaccepted(ptr_loc_state[1], A)
+ if isbounded(m) && !isacceptedLHA
# Add last value: the convention is that if the last transition is nothing,
# the trajectory is bounded
MarkovProcesses._finish_bounded_trajectory!(values, times, transitions, time_bound)
end
- if isabsorbing && !isacceptedLHA
+ if (isabsorbing || tn > time_bound) && !isacceptedLHA
next_state!(A, ptr_loc_state, values_state, ptr_time_state,
xn, time_bound, nothing, xn, p_sim, edge_candidates; verbose = verbose)
end
@@ -398,6 +398,7 @@ function simulate(pm::ParametricModel, p_prior::AbstractVector{Float64})
return simulate(pm.m; p = full_p)
end
+
"""
`volatile_simulate(pm::ParametricModel, p_prior::AbstractVector{Float64})
@@ -413,6 +414,25 @@ function volatile_simulate(pm::ParametricModel, p_prior::AbstractVector{Float64}
return volatile_simulate(pm.m; p = full_p, epsilon = epsilon)
end
+
+"""
+ `change_simulation_stop_criteria(m::ContinuousTimeModel, isabsorbing_func::Symbol)`
+
+ Change the simulation of the model `m` by adding a stop criteria based on the function named `isabsorbing_func::Symbol`.
+ isabsorbing_func must have the type signature `isabsorbing_func(p::Vector{Float64}, x::Vector{Int})` where `p` is the parameter vector of the model and `x` a state (not an observed state) of the model.
+"""
+function change_simulation_stop_criteria(m::ContinuousTimeModel, isabsorbing_func::Symbol)
+ model_name = Symbol(typeof(m))
+ isabs = getfield(Main, isabsorbing_func)
+ try
+ isabs(m.p, m.x0)
+ catch err
+ error("Something went wrong when trying to apply the criteria function")
+ end
+ @everywhere @eval $(generate_code_simulation(model_name, m.f!, isabsorbing_func; run_isabsorbing = true))
+end
+
+
"""
`distribute_mean_value_lha(sm::SynchronizedModel, sym_var::Symbol, nbr_stim::Int)`
diff --git a/tests/automata/euclidean_distance_on_bound.jl b/tests/automata/euclidean_distance_on_bound.jl
index dae30d13a66552c755f6798266f58eb683c67450..51ad64dd65dc2e4ace752b664c7151b658954698 100644
--- a/tests/automata/euclidean_distance_on_bound.jl
+++ b/tests/automata/euclidean_distance_on_bound.jl
@@ -6,36 +6,86 @@ import Distributions: Uniform
MAKE_SECOND_AUTOMATON_TESTS = false
load_model("SIR")
+load_automaton("euclidean_distance_automaton")
+if MAKE_SECOND_AUTOMATON_TESTS
+ load_automaton("euclidean_distance_automaton_2")
+end
+
tml_obs = 0:0.5:200
set_time_bound!(SIR, 200.0)
y_obs = vectorize(simulate(SIR), :I, tml_obs)
-
-load_automaton("euclidean_distance_automaton")
aut1 = create_euclidean_distance_automaton(SIR, tml_obs, y_obs, :I)
sync_SIR = SIR * aut1
-σ = simulate(sync_SIR)
-test = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
-if !test
- @show euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
- @show σ
+test_all = true
+nbr_sim = 100
+
+for i = 1:nbr_sim
+ let σ, test
+ σ = simulate(sync_SIR)
+ test = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
+ if !test
+ @show tml_obs
+ @show euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
+ @show σ
+ end
+ global test_all = test_all && test
+ end
end
if MAKE_SECOND_AUTOMATON_TESTS
- load_automaton("euclidean_distance_automaton_2")
aut2 = create_euclidean_distance_automaton_2(SIR, tml_obs, y_obs, :I)
sync_SIR = SIR * aut2
- σ = simulate(sync_SIR)
- test2 = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
- if !test2
- @show euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
- @show σ
+ for i = 1:nbr_sim
+ let σ, test2
+ σ = simulate(sync_SIR)
+ test2 = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
+ if !test2
+ @show euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
+ @show σ
+ end
+ global test_all = test_all && test2
+ end
+ end
+end
+
+tml_obs = 0:20.0:200
+set_time_bound!(SIR, 200.0)
+y_obs = vectorize(simulate(SIR), :I, tml_obs)
+aut1 = create_euclidean_distance_automaton(SIR, tml_obs, y_obs, :I)
+sync_SIR = SIR * aut1
+
+test_all = true
+nbr_sim = 100
+
+for i = 1:nbr_sim
+ let σ, test
+ σ = simulate(sync_SIR)
+ test = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
+ if !test
+ @show tml_obs
+ @show euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
+ @show σ
+ end
+ global test_all = test_all && test
end
-else
- test2 = true
end
-test_all = test && test2
+if MAKE_SECOND_AUTOMATON_TESTS
+ aut2 = create_euclidean_distance_automaton_2(SIR, tml_obs, y_obs, :I)
+ sync_SIR = SIR * aut2
+ for i = 1:nbr_sim
+ let σ, test2
+ σ = simulate(sync_SIR)
+ test2 = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
+ if !test2
+ @show euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
+ @show σ
+ end
+ global test_all = test_all && test2
+ end
+ end
+end
return test_all
diff --git a/tests/run_unit.jl b/tests/run_unit.jl
index deb8bb1bc1834dc9edc6820efbf8a2c1f118c685..21861b5ed7e81ef734c9b5b83ef2315850f29304 100644
--- a/tests/run_unit.jl
+++ b/tests/run_unit.jl
@@ -41,6 +41,7 @@ using Test
@test include("unit/simulate_sir.jl")
@test include("unit/simulate_sir_bounded.jl")
@test include("unit/simulate_er.jl")
+ @test include("unit/simulation_stop_criteria.jl")
@test include("unit/square_wave_oscillator.jl")
@test include("unit/vectorize.jl")
diff --git a/tests/unit/simulation_stop_criteria.jl b/tests/unit/simulation_stop_criteria.jl
new file mode 100644
index 0000000000000000000000000000000000000000..f20fb8fdedabb31ca6d2accb39c320dd03422b7f
--- /dev/null
+++ b/tests/unit/simulation_stop_criteria.jl
@@ -0,0 +1,62 @@
+
+using MarkovProcesses
+
+# Settings with two SIR models
+load_model("SIR")
+SIR2 = @network_model begin
+ R1: (S+I => 2I, ki*S*I)
+ R2: (I => R, kr*I)
+end "SIR2"
+set_x0!(SIR, [:S, :I, :R], [95, 5, 0])
+set_x0!(SIR2, [:S, :I, :R], [95, 5, 0])
+set_param!(SIR, [:ki, :kr], [0.012, 0.05])
+set_param!(SIR2, [:ki, :kr], [0.012, 0.05])
+# Settings ER model
+load_model("ER")
+set_x0!(ER, [:E, :S, :ES, :P], [100, 100, 0, 0])
+set_param!(ER, [:k1, :k2, :k3], [1.0, 1.0, 1.0])
+# Settings repressilator model
+load_model("repressilator")
+set_observed_var!(repressilator, [:mRNA1, :mRNA2, :mRNA3, :P1, :P2, :P3])
+set_x0!(repressilator, [:mRNA1, :mRNA2, :mRNA3], fill(0, 3))
+set_x0!(repressilator, [:P1, :P2, :P3], [5, 0, 15])
+set_param!(repressilator, :n, 2.0)
+set_param!(repressilator, [:α, :α0, :β, :n], [400.0, 0.0, 2.0, 2.0])
+set_time_bound!(repressilator, 200.0)
+
+test_all = true
+
+# First test : infinite simulation
+set_time_bound!(ER, Inf)
+set_time_bound!(SIR, Inf)
+set_time_bound!(SIR2, Inf)
+for i = 1:10
+ global test_all = test_all && simulate(ER).P[end] == 100
+ global test_all = test_all && simulate(SIR).I[end] == 0
+ global test_all = test_all && simulate(SIR2).I[end] == 0
+end
+
+# Second test: change the simulation stop criteria
+get_idx(m::ContinuousTimeModel, var::Symbol) = m.map_var_idx[var]
+idx_I_SIR, idx_I_SIR2, idx_P_ER, idx_P1_repressilator =
+get_idx(SIR, :I), get_idx(SIR2, :I), get_idx(ER, :P), get_idx(repressilator, :P1)
+
+@eval new_isabsorbing_SIR(p::Vector{Float64}, x::Vector{Int}) = x[$(idx_I_SIR)] == 2
+@eval new_isabsorbing_SIR2(p::Vector{Float64}, x::Vector{Int}) = x[$(idx_I_SIR2)] == 3
+@eval new_isabsorbing_ER(p::Vector{Float64}, x::Vector{Int}) = x[$(idx_P_ER)] == 50
+@eval new_isabsorbing_repressilator(p::Vector{Float64}, x::Vector{Int}) = x[$(idx_P1_repressilator)] == 100
+
+change_simulation_stop_criteria(SIR, :new_isabsorbing_SIR)
+change_simulation_stop_criteria(SIR2, :new_isabsorbing_SIR2)
+change_simulation_stop_criteria(ER, :new_isabsorbing_ER)
+change_simulation_stop_criteria(repressilator, :new_isabsorbing_repressilator)
+
+for i = 1:10
+ global test_all = test_all && simulate(SIR).I[end] == 2
+ global test_all = test_all && simulate(SIR2).I[end] == 3
+ global test_all = test_all && simulate(ER).P[end] == 50
+ global test_all = test_all && simulate(repressilator).P1[end] == 100
+end
+
+return test_all
+