Commit f98cd441 authored by Bentriou Mahmoud's avatar Bentriou Mahmoud

The new meta paradigm was extended to all automata and models

implemetend in the package.
A lot of files has been modified because some of the internal syntax has
changed. But the top level methods should still work.

Two tests about the plots of synchronized oscillatory simulations have
been added.

All tests passed.
parent dd619a95
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
# Creation of the automaton types
lha_name = :EuclideanDistanceAutomaton
edge_type = :EdgeEuclideanDistanceAutomaton
@everywhere @eval abstract type $(edge_type) <: Edge end
@everywhere @eval $(MarkovProcesses.generate_code_lha_type_def(lha_name, edge_type))
@everywhere @eval abstract type EdgeEuclideanDistanceAutomaton <: Edge end
@everywhere @eval $(MarkovProcesses.generate_code_lha_type_def(:EuclideanDistanceAutomaton, :EdgeEuclideanDistanceAutomaton))
function create_euclidean_distance_automaton(m::ContinuousTimeModel, timeline::AbstractVector{Float64}, observations::AbstractVector{Float64}, sym_obs::VariableModel)
# Requirements for the automaton
......@@ -47,15 +45,15 @@ function create_euclidean_distance_automaton(m::ContinuousTimeModel, timeline::A
Symbol("Edge_$(lha_name)_$(basename_func)_$(from_loc)$(to_loc)_$(edge_number)")
## check_constraints & update_state!
meta_elementary_func = quote
@everywhere @eval begin
# l0 loc
# l0 => l1
struct $(edge_name(:l0, :l1, 1)) <: $(edge_type) transitions::Union{Nothing,Vector{Symbol}} end
$(check_constraints)(edge::$(edge_name(:l0, :l1, 1)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) = true
$(update_state!)(edge::$(edge_name(:l0, :l1, 1)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
(setindex!(S_values, x[$(idx_obs_var)], $(to_idx(:n)));
setindex!(S_values, 0.0, $(to_idx(:d)));
setindex!(S_values, 1.0, $(to_idx(:idx)));
(S_values[$(to_idx(:n))] = x[$(idx_obs_var)];
S_values[$(to_idx(:d))] = 0.0;
S_values[$(to_idx(:idx))] = 1.0;
:l1)
# l1 loc
......@@ -69,15 +67,14 @@ function create_euclidean_distance_automaton(m::ContinuousTimeModel, timeline::A
$(update_state!)(edge::$(edge_name(:l1, :l1, 1)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
(y_obs = $(Tuple(observations));
y_obs_idx = y_obs[convert(Int, S_values[$(to_idx(:idx))])];
setindex!(S_values, S_values[$(to_idx(:d))] + (S_values[$(to_idx(:n))] - y_obs_idx)^2,
$(to_idx(:d)));
setindex!(S_values, S_values[$(to_idx(:idx))] + 1.0, $(to_idx(:idx)));
S_values[$(to_idx(:d))] = S_values[$(to_idx(:d))]+(S_values[$(to_idx(:n))]-y_obs_idx)^2;
S_values[$(to_idx(:idx))] = S_values[$(to_idx(:idx))]+1.0;
:l1)
struct $(edge_name(:l1, :l1, 2)) <: $(edge_type) transitions::Union{Nothing,Vector{Symbol}} end
$(check_constraints)(edge::$(edge_name(:l1, :l1, 2)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) = true
$(update_state!)(edge::$(edge_name(:l1, :l1, 2)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
(setindex!(S_values, x[$(idx_obs_var)], $(to_idx(:n)));
(S_values[$(to_idx(:n))] = x[$(idx_obs_var)];
:l1)
# l1 => l2
......@@ -85,10 +82,9 @@ function create_euclidean_distance_automaton(m::ContinuousTimeModel, timeline::A
$(check_constraints)(edge::$(edge_name(:l1, :l2, 1)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
S_values[$(to_idx(:idx))] >= ($nbr_observations + 1)
$(update_state!)(edge::$(edge_name(:l1, :l2, 1)), S_time::Float64, S_values::Vector{Float64}, x::Vector{Int}, p::Vector{Float64}) =
(setindex!(S_values, sqrt(S_values[$(to_idx(:d))]), $(to_idx(:d)));
(S_values[$(to_idx(:d))] = sqrt(S_values[$(to_idx(:d))]);
:l2)
end
@everywhere eval($(meta_elementary_func))
@eval begin
map_edges = Dict{Location,Dict{Location,Vector{$(edge_type)}}}()
......@@ -99,17 +95,17 @@ function create_euclidean_distance_automaton(m::ContinuousTimeModel, timeline::A
## Edges
# l0 loc
# l0 => l1
edge1 = getfield(Main, $(Meta.quot(edge_name(:l0, :l1, 1))))(nothing)
edge1 = $(edge_name(:l0, :l1, 1))(nothing)
map_edges[:l0][:l1] = [edge1]
# l1 loc
# l1 => l1
edge1 = getfield(Main, $(Meta.quot(edge_name(:l1, :l1, 1))))(nothing)
edge2 = getfield(Main, $(Meta.quot(edge_name(:l1, :l1, 2))))([:ALL])
edge1 = $(edge_name(:l1, :l1, 1))(nothing)
edge2 = $(edge_name(:l1, :l1, 2))([:ALL])
map_edges[:l1][:l1] = [edge1, edge2]
# l1 => l2
edge1 = getfield(Main, $(Meta.quot(edge_name(:l1, :l2, 1))))(nothing)
edge1 = $(edge_name(:l1, :l2, 1))(nothing)
map_edges[:l1][:l2] = [edge1]
end
......@@ -120,7 +116,7 @@ function create_euclidean_distance_automaton(m::ContinuousTimeModel, timeline::A
constants[Symbol("y_$(convert(Float64, i))")] = observations[i]
end
# Updating next_state!
# Updating types and simulation methods
@everywhere @eval $(MarkovProcesses.generate_code_synchronized_model_type_def(model_name, lha_name))
@everywhere @eval $(MarkovProcesses.generate_code_next_state(lha_name, edge_type, check_constraints, update_state!))
@everywhere @eval $(MarkovProcesses.generate_code_synchronized_simulation(model_name, lha_name, edge_type, m.f!, m.isabsorbing))
......
This diff is collapsed.
This diff is collapsed.
......@@ -25,7 +25,7 @@ export LHA, StateLHA, Edge, Location, VariableAutomaton
export +, -, δ, dist_lp, euclidean_distance
export get_obs_var, length_states, length_obs_var
export get_state_from_time, get_var_from_time, vectorize
export isbounded, times, transitions
export isbounded, states, times, transitions
export check_consistency, issteadystate, isaccepted
# LHA related methods
......
This diff is collapsed.
mutable struct BenchmarkModel <: ContinuousTimeModel
dim_state::Int # state space dim
dim_params::Int # parameter space dim
map_var_idx::Dict{VariableModel,Int} # maps variable str to index in the state space
_map_obs_var_idx::Dict{VariableModel,Int} # maps variable str to index in the observed state space
map_param_idx::Dict{ParameterModel,Int} # maps parameter str to index in the parameter space
transitions::Vector{Transition}
p::Vector{Float64}
x0::Vector{Int}
t0::Float64
f!::Function
g::Vector{VariableModel} # of dimension dim_obs_state
_g_idx::Vector{Int} # of dimension dim_obs_state
isabsorbing::Function
time_bound::Float64
buffer_size::Int
estim_min_states::Int
end
function BenchmarkModel(dim_state::Int, dim_params::Int, map_var_idx::Dict{VariableModel,Int},
map_param_idx::Dict{ParameterModel,Int}, transitions::Vector{<:Transition},
p::Vector{Float64}, x0::Vector{Int}, t0::Float64,
f!::Function, isabsorbing::Function;
g::Vector{VariableModel} = [var for var in keys(map_var_idx)], time_bound::Float64 = Inf,
buffer_size::Int = 10, estim_min_states::Int = 50)
dim_obs_state = length(g)
transitions = convert(Vector{Transition}, transitions)
_map_obs_var_idx = Dict()
_g_idx = Vector{Int}(undef, dim_obs_state)
for i = 1:dim_obs_state
_g_idx[i] = map_var_idx[g[i]] # = ( (g[i] = i-th obs var)::VariableModel => idx in state space )
_map_obs_var_idx[g[i]] = i
end
if length(methods(f!)) >= 2
@warn "You have possibly redefined a function Model.f! used in a previously instantiated model."
end
if length(methods(isabsorbing)) >= 2
@warn "You have possibly redefined a function Model.isabsorbing used in a previously instantiated model."
end
new_model = BenchmarkModel(dim_state, dim_params, map_var_idx, _map_obs_var_idx, map_param_idx, transitions,
p, x0, t0, f!, g, _g_idx, isabsorbing, time_bound, buffer_size, estim_min_states)
@assert check_consistency(new_model)
return new_model
end
struct OldTrajectory <: AbstractTrajectory
m::ContinuousTimeModel
m::BenchmarkModel
values::Matrix{Int}
times::Vector{Float64}
transitions::Vector{Union{Symbol,Nothing}}
......@@ -27,7 +72,7 @@ _get_value_row(σ::OldTrajectory, var::Symbol, idx::Int) =
# Model
function _simulate_col(m::ContinuousTimeModel)
function _simulate_col(m::BenchmarkModel)
# trajectory fields
full_values = zeros(m.dim_state, 0)
times = zeros(0)
......@@ -61,7 +106,7 @@ function _simulate_col(m::ContinuousTimeModel)
return OldTrajectory(m, values, times, transitions)
end
function _simulate_row(m::ContinuousTimeModel)
function _simulate_row(m::BenchmarkModel)
# trajectory fields
full_values = zeros(m.dim_state, 0)
times = zeros(0)
......@@ -96,7 +141,7 @@ function _simulate_row(m::ContinuousTimeModel)
end
function _simulate_col_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
function _simulate_col_buffer(m::BenchmarkModel; buffer_size::Int = 5)
# trajectory fields
full_values = zeros(m.dim_state, 0)
times = zeros(0)
......@@ -136,7 +181,7 @@ function _simulate_col_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
return OldTrajectory(m, values, times, transitions)
end
function _simulate_row_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
function _simulate_row_buffer(m::BenchmarkModel; buffer_size::Int = 5)
# trajectory fields
full_values = zeros(0, m.dim_state)
times = zeros(0)
......@@ -176,7 +221,7 @@ function _simulate_row_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
return OldTrajectory(m, values, times, transitions)
end
function _simulate_without_view(m::ContinuousTimeModel)
function _simulate_without_view(m::BenchmarkModel)
# trajectory fields
full_values = Matrix{Int}(undef, 1, m.dim_state)
full_values[1,:] = m.x0
......@@ -222,7 +267,7 @@ function _simulate_without_view(m::ContinuousTimeModel)
end
# With trajectory values in Matrix type
function _simulate_27d56(m::ContinuousTimeModel)
function _simulate_27d56(m::BenchmarkModel)
# trajectory fields
full_values = Matrix{Int}(undef, 1, m.dim_state)
full_values[1,:] = m.x0
......@@ -259,9 +304,9 @@ function _simulate_27d56(m::ContinuousTimeModel)
if isbounded(m)
#=
if times[end] > m.time_bound
full_values[end,:] = full_values[end-1,:]
times[end] = m.time_bound
transitions[end] = nothing
full_values[end,:] = full_values[end-1,:]
times[end] = m.time_bound
transitions[end] = nothing
else
end
=#
......@@ -274,7 +319,7 @@ function _simulate_27d56(m::ContinuousTimeModel)
end
function _simulate_d7458(m::ContinuousTimeModel)
function _simulate_d7458(m::BenchmarkModel)
# trajectory fields
full_values = Vector{Vector{Int}}(undef, m.dim_state)
for i = 1:m.dim_state full_values[i] = Int[m.x0[i]] end
......
......@@ -9,6 +9,7 @@ copy(S::StateLHA) = StateLHA(getfield(S, :A), getfield(S, :loc), copy(getfield(S
# From the variable automaton var symbol this function get the index in S.values
get_idx_var_automaton(S::StateLHA, var::VariableAutomaton) = getfield(getfield(S, :A), :map_var_automaton_idx)[var]
get_value(S::StateLHA, idx_var::Int) = getfield(S, :values)[idx_var]
get_state_lha_value(A::LHA, values::Vector{Float64}, var::VariableAutomaton) = values[A.map_var_automaton_idx[var]]
getindex(S::StateLHA, var::VariableAutomaton) = getindex(getfield(S, :values), get_idx_var_automaton(S, var))
setindex!(S::StateLHA, val::Float64, var::VariableAutomaton) = setindex!(getfield(S, :values), val, get_idx_var_automaton(S, var))
......@@ -166,9 +167,7 @@ function generate_code_next_state(lha_name::Symbol, edge_type::Symbol,
if verbose
println("##### Begin next_state!")
@show xnplus1, tnplus1, tr_nplus1
@show Sn
end
# In terms of values not reference, Snplus1 == Sn
# First, we check the asynchronous transitions
while true
turns += 1
......@@ -284,10 +283,6 @@ function generate_code_next_state(lha_name::Symbol, edge_type::Symbol,
end
=#
end
#=
setfield!(Snplus1, :loc, ptr_loc_state[1])
setfield!(Snplus1, :time, ptr_time_state[1])
=#
if verbose
println("##### End next_state!")
end
......@@ -303,18 +298,24 @@ function read_trajectory(A::LHA, σ::AbstractTrajectory; verbose = false)
p_sim = proba_model.p
l_t = times(σ)
l_tr = transitions(σ)
Sn = init_state(A_new, σ[1], l_t[1])
Snplus1 = copy(Sn)
edge_candidates = Vector{Edge}(undef, 2)
S = init_state(A, σ[1], l_t[1])
ptr_loc = [S.loc]
ptr_time = [S.time]
values = S.values
lha_name = Symbol(typeof(A))
edge_type = getfield(Main, Symbol("Edge$(lha_name)"))
edge_candidates = Vector{edge_type}(undef, 2)
if verbose println("Init: ") end
if verbose @show Sn end
for n in 2:length_states(σ)
next_state!(Snplus1, A_new, σ[n], l_t[n], l_tr[n], Sn, σ[n-1], p_sim, edge_candidates; verbose = verbose)
copyto!(Sn, Snplus1)
if Snplus1.loc in A_new.locations_final
if verbose @show S end
func_next_state! = getfield(Main, :next_state!) # It's defined after the compilation of the package
for k in 2:length_states(σ)
func_next_state!(A, ptr_loc, values, ptr_time, σ[k], l_t[k], l_tr[k], σ[k-1], p_sim, edge_candidates; verbose = verbose)
if ptr_loc[1] in A_new.locations_final
break
end
end
return Sn
setfield!(S, :loc, ptr_loc[1])
setfield!(S, :time, ptr_time[1])
return S
end
......@@ -385,10 +385,10 @@ The model in pm should be of type SynchronizedModel (`typeof(pm.m) <: Synchroniz
It returns `S::StateLHA`, not a trajectory.
"""
function volatile_simulate(pm::ParametricModel, p_prior::AbstractVector{Float64};
epsilon::Float64)
epsilon::Union{Nothing,Float64} = nothing)
@assert typeof(pm.m) <: SynchronizedModel
# ABC related automata
if typeof(pm.m.A) in <: EuclideanDistanceABCAutomaton
if @isdefined(EuclideanDistanceABCAutomaton) && typeof(pm.m.automaton) <: EuclideanDistanceABCAutomaton
nothing
end
full_p = copy(get_proba_model(pm).p)
......
......@@ -22,7 +22,7 @@ function plot(σ::AbstractTrajectory, vars::VariableModel...; plot_transitions::
end
# Plots
p = plot(title = "Trajectory of $(σ.m.name) model", palette = :lightrainbow, legend = :outertopright, background_color_legend=:transparent, dpi = 480)
p = plot(title = "Trajectory of $(Symbol(typeof(σ.m))) model", palette = :lightrainbow, legend = :outertopright, background_color_legend=:transparent, dpi = 480)
for var in to_plot
@assert var in get_obs_var(σ) "Variable $var is not observed."
plot!(p, times(σ), σ[var],
......@@ -55,24 +55,18 @@ end
function plot!(A::LHA; label::String = "")
label_region = (label == "") ? "Region LHA" : label
if A.name in ["F property", "G property"]
if haskey(A.constants, :x1) && haskey(A.constants, :x2) && haskey(A.constants, :t1) && haskey(A.constants, :t2)
x1, x2, t1, t2 = A.constants[:x1], A.constants[:x2], A.constants[:t1], A.constants[:t2]
plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5, label = label_region)
end
if (@isdefined(AutomatonF) && (typeof(A) <: AutomatonF)) || (@isdefined(AutomatonG) && (typeof(A) <: AutomatonG))
x1, x2, t1, t2 = A.constants[:x1], A.constants[:x2], A.constants[:t1], A.constants[:t2]
plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5, label = label_region)
end
if A.name in ["G and F property"]
if haskey(A.constants, :x1) && haskey(A.constants, :x2) && haskey(A.constants, :t1) && haskey(A.constants, :t2)
x1, x2, t1, t2 = A.constants[:x1], A.constants[:x2], A.constants[:t1], A.constants[:t2]
plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5, label = label_region)
end
if haskey(A.constants, :x3) && haskey(A.constants, :x4) && haskey(A.constants, :t3) && haskey(A.constants, :t4)
x3, x4, t3, t4 = A.constants[:x3], A.constants[:x4], A.constants[:t3], A.constants[:t4]
plot!(Shape([(t3,x3), (t3,x4), (t4,x4), (t4,x3), (t3,x3)]), opacity = 0.5, label = label_region)
end
if @isdefined(AutomatonGandF) && typeof(A) <: AutomatonGandF
x1, x2, t1, t2 = A.constants[:x1], A.constants[:x2], A.constants[:t1], A.constants[:t2]
plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5, label = label_region)
x3, x4, t3, t4 = A.constants[:x3], A.constants[:x4], A.constants[:t3], A.constants[:t4]
plot!(Shape([(t3,x3), (t3,x4), (t4,x4), (t4,x3), (t3,x3)]), opacity = 0.5, label = label_region)
end
label_region = (label == "") ? "L, H" : label
if A.name == "Period"
if @isdefined(PeriodAutomaton) && typeof(A) <: PeriodAutomaton
hline!([A.constants[:L], A.constants[:H]], label = label_region, linestyle = :dot)
end
end
......@@ -81,35 +75,38 @@ end
function plot_periodic_trajectory(A::LHA, σ::SynchronizedTrajectory, sym_obs::Symbol;
verbose = false, annot_size = 6, show_tp = false, filename::String = "")
@assert sym_obs in get_obs_var(σ) "Variable is not observed in the model"
@assert A.name in ["Period"]
@assert typeof(A) <: PeriodAutomaton "The automaton is not a period automaton"
@assert σ.m._g_idx == 1:σ.m.dim_state "All the variables must be observed. Call observe_all!(model) before."
p_sim = (σ.m).p
l_t = times(σ)
l_tr = transitions(σ)
Sn = init_state(A, σ[1], l_t[1])
Snplus1 = deepcopy(Sn)
S0 = init_state(A, σ[1], l_t[1])
ptr_loc = [S0.loc]
ptr_time = [S0.time]
values = S0.values
nbr_states = length_states(σ)
locations_trajectory = Vector{Location}(undef, nbr_states)
locations_trajectory[1] = Sn.loc
locations_trajectory[1] = S0.loc
idx_n = [1]
values_n = [Sn[:n]]
values_tp = [Sn[:tp]]
edge_candidates = Vector{Edge}(undef, 2)
values_n = [MarkovProcesses.get_state_lha_value(A, values, :n)]
values_tp = [MarkovProcesses.get_state_lha_value(A, values, :tp)]
edge_candidates = Vector{EdgePeriodAutomaton}(undef, 2)
if verbose println("Init: ") end
if verbose @show Sn end
for n in 2:nbr_states
next_state!(Snplus1, A, σ[n], l_t[n], l_tr[n], Sn, σ[n-1], p_sim, edge_candidates; verbose = verbose)
if Snplus1[:n] != values_n[end]
push!(idx_n, n)
push!(values_n, Snplus1[:n])
push!(values_tp, Sn[:tp])
for k in 2:nbr_states
tp_k = MarkovProcesses.get_state_lha_value(A, values, :tp)
next_state!(A, ptr_loc, values, ptr_time, σ[k], l_t[k], l_tr[k], σ[k-1], p_sim, edge_candidates; verbose = verbose)
n_kplus1 = MarkovProcesses.get_state_lha_value(A, values, :n)
if n_kplus1 != values_n[end]
push!(idx_n, k)
push!(values_n, n_kplus1)
push!(values_tp, tp_k)
end
copyto!(Sn, Snplus1)
locations_trajectory[n] = Sn.loc
if Snplus1.loc in A.locations_final
break
locations_trajectory[k] = ptr_loc[1]
if ptr_loc[1] in A.locations_final
break
end
end
p = plot(title = "Oscillatory trajectory of $(σ.m.name) model", palette = :lightrainbow,
p = plot(title = "Oscillatory trajectory of $(Symbol(typeof(σ.m))) model", palette = :lightrainbow,
background_color_legend=:transparent, dpi = 480, legend = :outertopright) #legendfontsize, legend
colors_loc = Dict(:l0 => :silver, :l0prime => :silver, :final => :black,
:low => :skyblue2, :mid => :orange, :high => :red)
......
......@@ -46,6 +46,6 @@ 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!,isabsorbing_ER_col; g=g)
ER_col = BenchmarkModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,ER_col_f!,isabsorbing_ER_col; g=g)
export ER_col
......@@ -45,6 +45,6 @@ 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!,isabsorbing_ER_col_buffer; g=g)
ER_col_buffer = BenchmarkModel(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
......@@ -45,6 +45,6 @@ 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!,isabsorbing_ER_row_buffer; g=g)
ER_row_buffer = BenchmarkModel(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
......@@ -44,6 +44,6 @@ end
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!,isabsorbing_SIR_col; g=g)
SIR_col = BenchmarkModel(d,k,dict_var,dict_p,l_tr,p,x0,t0,SIR_col_f!,isabsorbing_SIR_col; g=g)
export SIR_col
......@@ -44,6 +44,6 @@ end
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!,isabsorbing_SIR_col_buffer; g=g)
SIR_col_buffer = BenchmarkModel(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
......@@ -44,6 +44,6 @@ end
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!,isabsorbing_SIR_row_buffer; g=g)
SIR_row_buffer = BenchmarkModel(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
......@@ -6,7 +6,7 @@ doping_3way_oscillator = @network_model begin
R4: (DA + C => DA + A, rC*C*DA)
R5: (DB + A => DB + B, rA*A*DB)
R6: (DC + B => DC + C, rB*B*DC)
end "doping 3way oscillator"
end "Doping3wayOscillator"
set_x0!(doping_3way_oscillator, [:A,:B,:C,:DA,:DB,:DC], [333,333,333,10,10,10])
set_param!(doping_3way_oscillator, [:rA,:rB,:rC], [1.0, 1.0, 1.0])
......
......@@ -12,7 +12,7 @@ repressilator = @network_model begin
degr4: (P1 => 0, P1)
degr5: (P2 => 0, P2)
degr6: (P3 => 0, P3)
end "Repressilator pkg"
end "RepressilatorPkg"
set_observed_var!(repressilator, [:mRNA1, :mRNA2, :mRNA3, :P1, :P2, :P3])
set_x0!(repressilator, [:mRNA1, :mRNA2, :mRNA3], fill(0, 3))
......
......@@ -10,6 +10,7 @@ load_model("ER")
observe_all!(SIR)
observe_all!(ER)
MAKE_SECOND_AUTOMATON_TESTS = false
test_all = true
# SIR model
......@@ -30,16 +31,20 @@ for i = 1:nbr_sim
global sync_model = sync_SIR
break
end
sync_SIR = SIR * create_euclidean_distance_automaton_2(SIR, tml_obs, y_obs, :I)
σ = simulate(sync_SIR)
test2 = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
if !test2
@show test2, euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
global err = σ
global tml = tml_obs
global y = y_obs
global sync_model = sync_SIR
break
if MAKE_SECOND_AUTOMATON_TESTS
sync_SIR = SIR * create_euclidean_distance_automaton_2(SIR, tml_obs, y_obs, :I)
σ = simulate(sync_SIR)
test2 = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
if !test2
@show test2, euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
global err = σ
global tml = tml_obs
global y = y_obs
global sync_model = sync_SIR
break
end
else
test2 = true
end
global test_all = test_all && test && test2
end
......@@ -62,18 +67,21 @@ for i = 1:nbr_sim
global sync_model = sync_ER
break
end
sync_ER = ER * create_euclidean_distance_automaton_2(ER, tml_obs, y_obs, :P)
σ = simulate(sync_ER)
test2 = euclidean_distance(σ, :P, tml_obs, y_obs) == σ.state_lha_end[:d]
if !test2
@show test2, euclidean_distance(σ, :P, tml_obs, y_obs), σ.state_lha_end[:d]
global err = σ
global tml = tml_obs
global y = y_obs
global sync_model = sync_ER
break
if MAKE_SECOND_AUTOMATON_TESTS
sync_ER = ER * create_euclidean_distance_automaton_2(ER, tml_obs, y_obs, :P)
σ = simulate(sync_ER)
test2 = euclidean_distance(σ, :P, tml_obs, y_obs) == σ.state_lha_end[:d]
if !test2
@show test2, euclidean_distance(σ, :P, tml_obs, y_obs), σ.state_lha_end[:d]
global err = σ
global tml = tml_obs
global y = y_obs
global sync_model = sync_ER
break
end
else
test2 = true
end
#@show test, euclidean_distance(σ, tml_obs, y_obs, :P), σ.state_lha_end[:d]
global test_all = test_all && test && test2
end
end
......
......@@ -3,22 +3,28 @@ using MarkovProcesses
import LinearAlgebra: dot
import Distributions: Uniform
load_automaton("euclidean_distance_automaton")
load_automaton("euclidean_distance_automaton_2")
MAKE_SECOND_AUTOMATON_TESTS = false
load_model("SIR")
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]
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 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]