Commit 8041f629 authored by Bentriou Mahmoud's avatar Bentriou Mahmoud
Browse files

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.
parent 44feda9e
......@@ -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
......
......@@ -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
......
......@@ -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.Λ,
......
......@@ -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)
......
......@@ -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)
......
......@@ -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")
......
......@@ -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
......@@ -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
......@@ -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
......@@ -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
......@@ -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
......@@ -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
......@@ -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
......@@ -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
......@@ -5,8 +5,22 @@ 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
sync_SIR = SIR * A_F
σ, state_lha = simulate(sync_SIR)
return test_all
......@@ -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
......@@ -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
......
......@@ -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
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment