Add of a new feature: change_simulation_stop_criteria. It allows to

change the stop criteration of a ContinuousTimeModel simulation.
Minor fix  in ContinuousTimeModel for the euclidean distance automata
tests.
Add of a test for stop criteria simulation of ContinousTimeModel
(unit/simulation_stop_criteria).

All tests passed.

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!

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)`
 

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
 

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

tests/unit/simulation_stop_criteria.jl

+
+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
+