From 8041f629c1cc0c4333ddfcfb6938e905f5b225a6 Mon Sep 17 00:00:00 2001 From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr> Date: Sat, 21 Nov 2020 20:01:35 +0100 Subject: [PATCH] Rename of functions is_X => isX because of Julia Base name functions. Add of a new type of Trajectory: SynchronizedTrajectory which is created by the simulation of a synchronized model with an LHA. Simulate function for SynchronizedModel executes and is tested by two simple models. Next step : Cosmos based tests. --- core/MarkovProcesses.jl | 12 +-- core/_tests_simulate.jl | 52 ++++++------- core/common.jl | 20 +++-- core/lha.jl | 2 +- core/model.jl | 77 ++++++++++---------- core/trajectory.jl | 9 ++- models/ER.jl | 4 +- models/SIR.jl | 4 +- models/_bench_perf_test/ER_col.jl | 4 +- models/_bench_perf_test/ER_col_buffer.jl | 4 +- models/_bench_perf_test/ER_row_buffer.jl | 4 +- models/_bench_perf_test/SIR_col.jl | 4 +- models/_bench_perf_test/SIR_col_buffer.jl | 4 +- models/_bench_perf_test/SIR_row_buffer.jl | 4 +- tests/automata/automaton_F.jl | 12 --- tests/automata/sync_simulate_last_state_F.jl | 26 +++++++ tests/run_automata.jl | 1 + tests/unit/is_always_bounded_sir.jl | 2 +- tests/unit/side_effects_1.jl | 2 +- 19 files changed, 137 insertions(+), 110 deletions(-) delete mode 100644 tests/automata/automaton_F.jl create mode 100644 tests/automata/sync_simulate_last_state_F.jl diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl index 08c2b4a..4ae10bb 100644 --- a/core/MarkovProcesses.jl +++ b/core/MarkovProcesses.jl @@ -6,23 +6,23 @@ import Base: copy, getfield, getindex, lastindex, setindex!, getproperty, setpro import StaticArrays: SVector # Common types and constructors -export Observations, AbstractTrajectory, Trajectory +export Observations, AbstractTrajectory, Trajectory, SynchronizedTrajectory +export Model, ContinuousTimeModel, SynchronizedModel export LHA, StateLHA, Edge -export Model, ContinuousTimeModel, DiscreteTimeModel # Trajectory related methods export +, -, δ, dist_lp export get_obs_var, length_states, length_obs_var, get_state_from_time -export is_bounded, times, transitions -export check_consistency, is_steadystate +export isbounded, times, transitions +export check_consistency, issteadystate, isaccepted # LHA related methods export init_state, next_state!, read_trajectory -export load_automaton, get_index, get_value, length_var +export load_automaton, get_index, get_value, length_var, isaccepted # Model related methods export simulate, set_param!, get_param, set_observed_var! -export is_bounded, check_consistency +export isbounded, isaccepted, check_consistency export load_model, get_module_path # Utils diff --git a/core/_tests_simulate.jl b/core/_tests_simulate.jl index 0fc7c4a..632b3ff 100644 --- a/core/_tests_simulate.jl +++ b/core/_tests_simulate.jl @@ -33,18 +33,18 @@ function _simulate_col(m::ContinuousTimeModel) # at time n+1 xnplus1 = zeros(Int, m.d) tnplus1 = zeros(Float64, 1) - is_absorbing = (m.is_absorbing(m.p,xn))::Bool - while !is_absorbing && (tn <= m.time_bound) + isabsorbing = (m.isabsorbing(m.p,xn))::Bool + while !isabsorbing && (tn <= m.time_bound) m.f!(xnplus1, tnplus1, tr, xn, tn, m.p) full_values = hcat(full_values, xnplus1) push!(times, tnplus1[1]) push!(transitions, tr[1]) xn, tn = xnplus1, tnplus1[1] n += 1 - is_absorbing = m.is_absorbing(m.p,xn)::Bool + isabsorbing = m.isabsorbing(m.p,xn)::Bool end values = @view full_values[m._g_idx,:] - if is_bounded(m) + if isbounded(m) if times[end] > m.time_bound values[:, end] = values[:,end-1] times[end] = m.time_bound @@ -67,18 +67,18 @@ function _simulate_row(m::ContinuousTimeModel) # at time n+1 xnplus1 = zeros(Int, m.d) tnplus1 = zeros(Float64, 1) - is_absorbing = (m.is_absorbing(m.p,xn))::Bool - while !is_absorbing && (tn <= m.time_bound) + isabsorbing = (m.isabsorbing(m.p,xn))::Bool + while !isabsorbing && (tn <= m.time_bound) m.f!(xnplus1, tnplus1, tr, xn, tn, m.p) full_values = vcat(full_values, xnplus1) push!(times, tnplus1[1]) push!(transitions, tr[1]) xn, tn = xnplus1, tnplus1[1] n += 1 - is_absorbing = m.is_absorbing(m.p,xn)::Bool + isabsorbing = m.isabsorbing(m.p,xn)::Bool end values = @view full_values[m._g_idx,:] - if is_bounded(m) + if isbounded(m) if times[end] > m.time_bound values[:, end] = values[:,end-1] times[end] = m.time_bound @@ -102,24 +102,24 @@ function _simulate_col_buffer(m::ContinuousTimeModel; buffer_size::Int = 5) mat_x = zeros(Int, m.d, buffer_size) l_t = zeros(Float64, buffer_size) l_tr = Vector{String}(undef, buffer_size) - is_absorbing = m.is_absorbing(m.p,xn)::Bool - while !is_absorbing && (tn <= m.time_bound) + isabsorbing = m.isabsorbing(m.p,xn)::Bool + while !isabsorbing && (tn <= m.time_bound) i = 0 - while i < buffer_size && !is_absorbing && (tn <= m.time_bound) + while i < buffer_size && !isabsorbing && (tn <= m.time_bound) i += 1 m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p) xn = @view mat_x[:,i] tn = l_t[i] - is_absorbing = m.is_absorbing(m.p,xn)::Bool + isabsorbing = m.isabsorbing(m.p,xn)::Bool end full_values = hcat(full_values, @view mat_x[:,1:i]) append!(times, @view l_t[1:i]) append!(transitions, @view l_tr[1:i]) n += i - is_absorbing = m.is_absorbing(m.p,xn)::Bool + isabsorbing = m.isabsorbing(m.p,xn)::Bool end values = @view full_values[m._g_idx,:] - if is_bounded(m) + if isbounded(m) if times[end] > m.time_bound values[:, end] = values[:,end-1] times[end] = m.time_bound @@ -142,24 +142,24 @@ function _simulate_row_buffer(m::ContinuousTimeModel; buffer_size::Int = 5) mat_x = zeros(Int, buffer_size, m.d) l_t = zeros(Float64, buffer_size) l_tr = Vector{String}(undef, buffer_size) - is_absorbing = m.is_absorbing(m.p,xn)::Bool - while !is_absorbing && (tn <= m.time_bound) + isabsorbing = m.isabsorbing(m.p,xn)::Bool + while !isabsorbing && (tn <= m.time_bound) i = 0 - while i < buffer_size && !is_absorbing && (tn <= m.time_bound) + while i < buffer_size && !isabsorbing && (tn <= m.time_bound) i += 1 m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p) xn = @view mat_x[i,:] tn = l_t[i] - is_absorbing = m.is_absorbing(m.p,xn)::Bool + isabsorbing = m.isabsorbing(m.p,xn)::Bool end full_values = vcat(full_values, @view mat_x[1:i,:]) append!(times, @view l_t[1:i]) append!(transitions, @view l_tr[1:i]) n += i - is_absorbing = m.is_absorbing(m.p,xn)::Bool + isabsorbing = m.isabsorbing(m.p,xn)::Bool end values = @view full_values[:,m._g_idx] - if is_bounded(m) + if isbounded(m) if times[end] > m.time_bound values[end,:] = values[end-1,:] times[end] = m.time_bound @@ -183,23 +183,23 @@ function _simulate_without_view(m::ContinuousTimeModel) mat_x = zeros(Int, m.buffer_size, m.d) l_t = zeros(Float64, m.buffer_size) l_tr = Vector{String}(undef, m.buffer_size) - is_absorbing = m.is_absorbing(m.p,xn)::Bool - while !is_absorbing && (tn <= m.time_bound) + isabsorbing = m.isabsorbing(m.p,xn)::Bool + while !isabsorbing && (tn <= m.time_bound) i = 0 - while i < m.buffer_size && !is_absorbing && (tn <= m.time_bound) + while i < m.buffer_size && !isabsorbing && (tn <= m.time_bound) i += 1 m.f!(mat_x, l_t, l_tr, i, @view(xn[:]), tn, m.p) xn = mat_x[i,:] tn = l_t[i] - is_absorbing = m.is_absorbing(m.p,@view(xn[:]))::Bool + isabsorbing = m.isabsorbing(m.p,@view(xn[:]))::Bool end full_values = vcat(full_values, mat_x[1:i,:]) append!(times, l_t[1:i]) append!(transitions, l_tr[1:i]) n += i - is_absorbing = m.is_absorbing(m.p,@view(xn[:]))::Bool + isabsorbing = m.isabsorbing(m.p,@view(xn[:]))::Bool end - if is_bounded(m) + if isbounded(m) if times[end] > m.time_bound full_values[end,:] = full_values[end-1,:] times[end] = m.time_bound diff --git a/core/common.jl b/core/common.jl index 2e9decc..bbe0914 100644 --- a/core/common.jl +++ b/core/common.jl @@ -19,7 +19,7 @@ mutable struct ContinuousTimeModel <: Model f!::Function g::Vector{String} # of dimension dobs _g_idx::Vector{Int} # of dimension dobs - is_absorbing::Function + isabsorbing::Function time_bound::Float64 buffer_size::Int end @@ -57,15 +57,23 @@ mutable struct StateLHA time::Float64 end -mutable struct SynchronizedModel +mutable struct SynchronizedModel <: Model m::ContinuousTimeModel automaton::LHA end +struct SynchronizedTrajectory <: AbstractTrajectory + S::StateLHA + m::SynchronizedModel + values::Matrix{Int} + times::Vector{Float64} + transitions::Vector{Transition} +end + # Constructors function ContinuousTimeModel(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_transitions::Vector{String}, p::Vector{Float64}, x0::Vector{Int}, t0::Float64, - f!::Function, is_absorbing::Function; + f!::Function, isabsorbing::Function; g::Vector{String} = keys(map_var_idx), time_bound::Float64 = Inf, buffer_size::Int = 10) dobs = length(g) _map_obs_var_idx = Dict() @@ -78,11 +86,11 @@ function ContinuousTimeModel(d::Int, k::Int, map_var_idx::Dict, map_param_idx::D if length(methods(f!)) >= 2 @warn "You have possibly redefined a function Model.f! used in a previously instantiated model." end - if length(methods(is_absorbing)) >= 2 - @warn "You have possibly redefined a function Model.is_absorbing used in a previously instantiated model." + if length(methods(isabsorbing)) >= 2 + @warn "You have possibly redefined a function Model.isabsorbing used in a previously instantiated model." end - return ContinuousTimeModel(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_transitions, p, x0, t0, f!, g, _g_idx, is_absorbing, time_bound, buffer_size) + return ContinuousTimeModel(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_transitions, p, x0, t0, f!, g, _g_idx, isabsorbing, time_bound, buffer_size) end LHA(A::LHA, map_var::Dict{String,Int}) = LHA(A.l_transitions, A.l_loc, A.Λ, diff --git a/core/lha.jl b/core/lha.jl index 3221688..c0288d6 100644 --- a/core/lha.jl +++ b/core/lha.jl @@ -19,7 +19,7 @@ function Base.show(io::IO, S::StateLHA) end end -isaccepted(S::StateLHA) = S.loc in (S.A).l_loc_final +isaccepted(S::StateLHA) = (S.loc in (S.A).l_loc_final) # Methods for synchronize / read the trajectory function init_state(A::LHA, x0::SubArray{Int,1}, t0::Float64) diff --git a/core/model.jl b/core/model.jl index 6824680..03099d0 100644 --- a/core/model.jl +++ b/core/model.jl @@ -16,32 +16,38 @@ function simulate(m::ContinuousTimeModel) mat_x = zeros(Int, m.buffer_size, m.d) l_t = zeros(Float64, m.buffer_size) l_tr = Vector{String}(undef, m.buffer_size) - is_absorbing::Bool = m.is_absorbing(m.p,xn) - while !is_absorbing && (tn <= m.time_bound) + isabsorbing::Bool = m.isabsorbing(m.p,xn) + while !isabsorbing && tn <= m.time_bound i = 0 - while i < m.buffer_size && !is_absorbing && (tn <= m.time_bound) + while i < m.buffer_size && !isabsorbing && tn <= m.time_bound i += 1 m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p) - xn = view(mat_x, i, :) tn = l_t[i] - is_absorbing = m.is_absorbing(m.p,xn) + if tn > m.time_bound + i -= 1 # 0 is an ok value, 1:0 is allowed + break + end + xn = view(mat_x, i, :) + isabsorbing = m.isabsorbing(m.p,xn) end full_values = vcat(full_values, view(mat_x, 1:i, :)) append!(times, view(l_t, 1:i)) append!(transitions, view(l_tr, 1:i)) n += i - is_absorbing = m.is_absorbing(m.p,xn) + #isabsorbing = m.isabsorbing(m.p,xn) end - if is_bounded(m) + if isbounded(m) + #= if times[end] > m.time_bound full_values[end,:] = full_values[end-1,:] times[end] = m.time_bound transitions[end] = nothing else - full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d)) - push!(times, m.time_bound) - push!(transitions, nothing) end + =# + full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d)) + push!(times, m.time_bound) + push!(transitions, nothing) end values = view(full_values, :, m._g_idx) return Trajectory(m, values, times, transitions) @@ -56,7 +62,7 @@ function simulate(product::SynchronizedModel) 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, reshaped_x0, t0) + S0 = init_state(A, reshaped_x0, m.t0) # values at time n n = 0 xn = reshaped_x0 @@ -66,47 +72,44 @@ function simulate(product::SynchronizedModel) mat_x = zeros(Int, m.buffer_size, m.d) l_t = zeros(Float64, m.buffer_size) l_tr = Vector{String}(undef, m.buffer_size) - is_absorbing::Bool = m.is_absorbing(m.p,xn) + isabsorbing::Bool = m.isabsorbing(m.p,xn) + isacceptedLHA::Bool = isaccepted(Sn) Snplus1 = copy(Sn) - while !is_absorbing && (tn < m.time_bound) + Sn_dump = Sn + while !isabsorbing && tn <= m.time_bound && !isacceptedLHA i = 0 - while i < m.buffer_size && !is_absorbing && (tn < m.time_bound) + 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) - xn = view(mat_x, i, :) tn = l_t[i] - tr_n = l_tr[n] - A.next_state!(Snplus1, A, xn, tn, tr_n, Sn) + 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 - is_absorbing = m.is_absorbing(m.p,xn) + isabsorbing = m.isabsorbing(m.p,xn) + isacceptedLHA = isaccepted(Snplus1) end full_values = vcat(full_values, view(mat_x, 1:i, :)) append!(times, view(l_t, 1:i)) append!(transitions, view(l_tr, 1:i)) n += i - is_absorbing = m.is_absorbing(m.p,xn) end - if is_bounded(m) - if times[end] > m.time_bound - full_values[end,:] = full_values[end-1,:] - times[end] = m.time_bound - transitions[end] = nothing - else - full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d)) - push!(times, m.time_bound) - push!(transitions, nothing) - end + # When the trajectory is accepted, we should not add an end value + if isbounded(m) && !isaccepted(Sn) + @assert times[end] < m.time_bound + #full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d)) + #push!(times, m.time_bound) + #push!(transitions, nothing) end values = view(full_values, :, m._g_idx) - return Trajectory(m, values, times, transitions) + return SynchronizedTrajectory(Snplus1, product, values, times, transitions) end function simulate(m::ContinuousTimeModel, n::Int) - obs = ContinuousObservations(undef, n) - for i = 1:n - obs[i] = simulate(m) - end - return obs end function set_observed_var!(m::Model,g::Vector{String}) @@ -122,7 +125,7 @@ function set_observed_var!(m::Model,g::Vector{String}) m._map_obs_var_idx = _map_obs_var_idx end -is_bounded(m::ContinuousTimeModel) = m.time_bound < Inf +isbounded(m::ContinuousTimeModel) = m.time_bound < Inf function check_consistency(m::ContinuousTimeModel) @assert m.d == length(m.map_var_idx) @assert m.k == length(m.map_param_idx) @@ -130,7 +133,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.is_absorbing(m.p, view(reshape(m.x0, 1, m.d), 1, :))) == Bool + @assert typeof(m.isabsorbing(m.p, view(reshape(m.x0, 1, m.d), 1, :))) == Bool return true end set_param!(m::ContinuousTimeModel, p::Vector{Float64}) = (m.p = p) diff --git a/core/trajectory.jl b/core/trajectory.jl index 32afb29..cb8901b 100644 --- a/core/trajectory.jl +++ b/core/trajectory.jl @@ -57,7 +57,7 @@ Requires `get_obs_var(σ1) == get_obs_var(σ2)`, it is verified if they are simu """ function dist_lp(σ1::AbstractTrajectory, σ2::AbstractTrajectory, var::String; verbose::Bool = false, p::Int = 1) - if !is_bounded(σ1) || !is_bounded(σ2) + if !isbounded(σ1) || !isbounded(σ2) @warn "Lp distance computed on unbounded trajectories. Result should be wrong" end return dist_lp(σ1[var], times(σ1), σ2[var], times(σ2); verbose = false, p = p) @@ -154,13 +154,14 @@ function check_consistency(σ::AbstractTrajectory) @assert length_obs_var(σ) == length(σ.m.g) == size(σ.values)[2] return true end -is_steadystate(σ::AbstractTrajectory) = (σ.m).is_absorbing((σ.m).p, σ[end]) +issteadystate(σ::AbstractTrajectory) = (σ.m).isabsorbing((σ.m).p, σ[end]) # Properties of the trajectory length_states(σ::AbstractTrajectory) = length(σ.times) length_obs_var(σ::AbstractTrajectory) = size(σ.values)[2] get_obs_var(σ::AbstractTrajectory) = (σ.m).g -is_bounded(σ::AbstractTrajectory) = σ.transitions[end] == nothing +isbounded(σ::AbstractTrajectory) = σ.transitions[end] == nothing +isaccepted(σ::SynchronizedTrajectory) = isaccepted(σ.S) # Access to trajectory values get_var_values(σ::AbstractTrajectory, var::String) = view(σ.values, :, (σ.m)._map_obs_var_idx[var]) @@ -172,7 +173,7 @@ function get_state_from_time(σ::AbstractTrajectory, t::Float64) l_t = times(σ) if t == l_t[end] return σ[end] end if t > l_t[end] - if !is_bounded(σ) + if !isbounded(σ) return σ[end] else error("This trajectory is bounded and you're accessing out of the bounds") diff --git a/models/ER.jl b/models/ER.jl index b427f08..deace15 100644 --- a/models/ER.jl +++ b/models/ER.jl @@ -41,10 +41,10 @@ function ER_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, i l_t[idx] = tn + tau l_tr[idx] = "R$(reaction)" end -is_absorbing_ER(p::Vector{Float64},xn::SubArray{Int,1}) = +isabsorbing_ER(p::Vector{Float64},xn::SubArray{Int,1}) = (p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0 g_ER = ["P"] -ER = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_ER,p_ER,x0_ER,t0_ER,ER_f!,is_absorbing_ER; g=g_ER) +ER = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_ER,p_ER,x0_ER,t0_ER,ER_f!,isabsorbing_ER; g=g_ER) export ER diff --git a/models/SIR.jl b/models/SIR.jl index d3fff55..37990bc 100644 --- a/models/SIR.jl +++ b/models/SIR.jl @@ -41,9 +41,9 @@ function SIR_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector{String}, l_t[idx] = tn + tau l_tr[idx] = "R$(reaction)" end -is_absorbing_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::SubArray{Int,1}) = (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!,is_absorbing_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) export SIR diff --git a/models/_bench_perf_test/ER_col.jl b/models/_bench_perf_test/ER_col.jl index 2546fd4..e548336 100644 --- a/models/_bench_perf_test/ER_col.jl +++ b/models/_bench_perf_test/ER_col.jl @@ -42,10 +42,10 @@ function ER_col_f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{St tnplus1[1] = tn + tau tr[1] = "R$(reaction)" end -is_absorbing_ER_col(p::Vector{Float64},xn::AbstractVector{Int}) = +isabsorbing_ER_col(p::Vector{Float64},xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0 g = ["P"] -ER_col = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_col_f!,is_absorbing_ER_col; g=g) +ER_col = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_col_f!,isabsorbing_ER_col; g=g) export ER_col diff --git a/models/_bench_perf_test/ER_col_buffer.jl b/models/_bench_perf_test/ER_col_buffer.jl index c7315c4..4261271 100644 --- a/models/_bench_perf_test/ER_col_buffer.jl +++ b/models/_bench_perf_test/ER_col_buffer.jl @@ -41,10 +41,10 @@ function ER_col_buffer_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector l_t[idx] = tn + tau l_tr[idx] = "R$(reaction)" end -is_absorbing_ER_col_buffer(p::Vector{Float64},xn::AbstractVector{Int}) = +isabsorbing_ER_col_buffer(p::Vector{Float64},xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0 g = ["P"] -ER_col_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_col_buffer_f!,is_absorbing_ER_col_buffer; g=g) +ER_col_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_col_buffer_f!,isabsorbing_ER_col_buffer; g=g) export ER_col_buffer diff --git a/models/_bench_perf_test/ER_row_buffer.jl b/models/_bench_perf_test/ER_row_buffer.jl index 5c1489d..0fb918f 100644 --- a/models/_bench_perf_test/ER_row_buffer.jl +++ b/models/_bench_perf_test/ER_row_buffer.jl @@ -41,10 +41,10 @@ function ER_row_buffer_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vector l_t[idx] = tn + tau l_tr[idx] = "R$(reaction)" end -is_absorbing_ER_row_buffer(p::Vector{Float64},xn::AbstractVector{Int}) = +isabsorbing_ER_row_buffer(p::Vector{Float64},xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0 g = ["P"] -ER_row_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_row_buffer_f!,is_absorbing_ER_row_buffer; g=g) +ER_row_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_row_buffer_f!,isabsorbing_ER_row_buffer; g=g) export ER_row_buffer diff --git a/models/_bench_perf_test/SIR_col.jl b/models/_bench_perf_test/SIR_col.jl index 3052b5e..381a487 100644 --- a/models/_bench_perf_test/SIR_col.jl +++ b/models/_bench_perf_test/SIR_col.jl @@ -41,9 +41,9 @@ function SIR_col_f!(xnplus1::Vector{Int}, tnplus1::Vector{Float64}, tr::Vector{S tnplus1[1] = tn + tau tr[1] = "R$(reaction)" end -is_absorbing_SIR_col(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0 +isabsorbing_SIR_col(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0 g = ["I"] -SIR_col = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_f!,is_absorbing_SIR_col; g=g) +SIR_col = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_f!,isabsorbing_SIR_col; g=g) export SIR_col diff --git a/models/_bench_perf_test/SIR_col_buffer.jl b/models/_bench_perf_test/SIR_col_buffer.jl index ee70393..771149d 100644 --- a/models/_bench_perf_test/SIR_col_buffer.jl +++ b/models/_bench_perf_test/SIR_col_buffer.jl @@ -41,9 +41,9 @@ function SIR_col_buffer_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vecto l_t[idx] = tn + tau l_tr[idx] = "R$(reaction)" end -is_absorbing_SIR_col_buffer(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0 +isabsorbing_SIR_col_buffer(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0 g = ["I"] -SIR_col_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_buffer_f!,is_absorbing_SIR_col_buffer; g=g) +SIR_col_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_buffer_f!,isabsorbing_SIR_col_buffer; g=g) export SIR_col_buffer diff --git a/models/_bench_perf_test/SIR_row_buffer.jl b/models/_bench_perf_test/SIR_row_buffer.jl index 4eaed57..ee241e6 100644 --- a/models/_bench_perf_test/SIR_row_buffer.jl +++ b/models/_bench_perf_test/SIR_row_buffer.jl @@ -41,9 +41,9 @@ function SIR_row_buffer_f!(mat_x::Matrix{Int}, l_t::Vector{Float64}, l_tr::Vecto l_t[idx] = tn + tau l_tr[idx] = "R$(reaction)" end -is_absorbing_SIR_row_buffer(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0 +isabsorbing_SIR_row_buffer(p::Vector{Float64}, xn::AbstractVector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0 g = ["I"] -SIR_row_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_row_buffer_f!,is_absorbing_SIR_row_buffer; g=g) +SIR_row_buffer = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_row_buffer_f!,isabsorbing_SIR_row_buffer; g=g) export SIR_row_buffer diff --git a/tests/automata/automaton_F.jl b/tests/automata/automaton_F.jl deleted file mode 100644 index 145da95..0000000 --- a/tests/automata/automaton_F.jl +++ /dev/null @@ -1,12 +0,0 @@ - -using MarkovProcesses - -load_model("SIR") -load_automaton("automaton_F") -SIR.time_bound = 120.0 -x1, x2, t1, t2 = 0.0, Inf, 100.0, 120.0 -A_F = create_automaton_F(SIR, x1, x2, t1, t2, "I") # <: LHA - -sync_SIR = SIR * A_F -σ, state_lha = simulate(sync_SIR) - diff --git a/tests/automata/sync_simulate_last_state_F.jl b/tests/automata/sync_simulate_last_state_F.jl new file mode 100644 index 0000000..7f65ae6 --- /dev/null +++ b/tests/automata/sync_simulate_last_state_F.jl @@ -0,0 +1,26 @@ + +using MarkovProcesses + +load_model("SIR") +load_automaton("automaton_F") +SIR.time_bound = 120.0 +x1, x2, t1, t2 = 0.0, Inf, 100.0, 120.0 + +A_F = create_automaton_F(SIR, x1, x2, t1, t2, "I") # <: LHA +sync_SIR = A_F * SIR + +function test_last_state(A::LHA, m::SynchronizedModel) + σ = simulate(m) + test = (get_state_from_time(σ, (σ.S).time)[1] == (σ.S)["n"]) && ((σ.S)["d"] == 0) + return test +end + +test_all = true +nbr_sim = 10000 +for i = 1:nbr_sim + test = test_last_state(A_F, sync_SIR) + global test_all = test_all && test +end + +return test_all + diff --git a/tests/run_automata.jl b/tests/run_automata.jl index abec4d0..d0b38fd 100644 --- a/tests/run_automata.jl +++ b/tests/run_automata.jl @@ -3,5 +3,6 @@ using Test @testset "Automata tests" begin @test include("automata/read_trajectory_last_state_F.jl") + @test include("automata/sync_simulate_last_state_F.jl") end diff --git a/tests/unit/is_always_bounded_sir.jl b/tests/unit/is_always_bounded_sir.jl index 99b4661..7b478e0 100644 --- a/tests/unit/is_always_bounded_sir.jl +++ b/tests/unit/is_always_bounded_sir.jl @@ -6,7 +6,7 @@ load_model("SIR") SIR.time_bound = 120.0 for i = 1:1000 σ = simulate(SIR) - global test = test && is_bounded(σ) + global test = test && isbounded(σ) end return test diff --git a/tests/unit/side_effects_1.jl b/tests/unit/side_effects_1.jl index 3392843..cd45dfb 100644 --- a/tests/unit/side_effects_1.jl +++ b/tests/unit/side_effects_1.jl @@ -10,7 +10,7 @@ test = SIR.map_var_idx !== ER.map_var_idx && SIR.g !== ER.g && SIR._g_idx !== ER._g_idx && SIR.f! != ER.f! && - SIR.is_absorbing != ER.is_absorbing + SIR.isabsorbing != ER.isabsorbing return test -- GitLab