Commit 95c83834 authored by Bentriou Mahmoud's avatar Bentriou Mahmoud

small fix in abc with reject automaton. add of more tests and benchmark

parent 5b3cbc2b
......@@ -19,9 +19,11 @@ end
map_var_automaton_idx::Dict{VariableAutomaton,Int} # nvar keys : str_var => idx in values
flow::Dict{Location,Vector{Float64}} # output of length nvar
map_edges::Dict{Location, Dict{Location,Vector{EdgeABCEuclideanDistanceAutomaton}}}
constants::Dict{Symbol,Float64}
constants::Dict
map_var_model_idx::Dict{VariableModel,Int} # of dim d (of a model)
ϵ::Float64
timeline::Vector{Float64}
observations::Vector{Float64}
end
# A push! method implementend by myself because of preallocation of edge_candidates
......@@ -336,18 +338,20 @@ function create_abc_euclidean_distance_automaton(m::ContinuousTimeModel, timelin
end
## Constants
constants = Dict{Symbol,Float64}(:nbr_obs => nbr_observations)
constants = Dict(:nbr_obs => nbr_observations, :species => sym_obs)
#=
for i = 1:nbr_observations
constants[Symbol("tml_$(convert(Float64, i))")] = timeline[i]
constants[Symbol("y_$(convert(Float64, i))")] = observations[i]
end
=#
# Updating types and simulation methods
@everywhere @eval $(MarkovProcesses.generate_code_synchronized_model_type_def(model_name, lha_name))
@everywhere @eval $(MarkovProcesses.generate_code_synchronized_simulation(model_name, lha_name, edge_type, m.f!, m.isabsorbing))
A = ABCEuclideanDistanceAutomaton(m.transitions, locations, Λ_F, locations_init, locations_final,
map_var_automaton_idx, flow, map_edges, constants, m.map_var_idx, Inf)
map_var_automaton_idx, flow, map_edges, constants, m.map_var_idx, Inf, timeline, observations)
return A
end
......
......@@ -8,6 +8,7 @@ import Distributions: Uniform
load_automaton("euclidean_distance_automaton")
load_automaton("euclidean_distance_automaton_2")
load_automaton("abc_euclidean_distance_automaton")
load_model("repressilator")
tb = 210.0
tml_obs = 0:10.0:210.0
......@@ -36,6 +37,23 @@ b_vol_sim_aut1 = @benchmark (σ = volatile_simulate($(sync1)))
@btime (σ = volatile_simulate($(sync1)))
@show minimum(b_vol_sim_aut1), mean(b_vol_sim_aut1), maximum(b_vol_sim_aut1)
println("ABC reject automaton with 1 loc")
aut1_abc = create_abc_euclidean_distance_automaton(repressilator, tml_obs, y_obs, :P1)
aut1_abc.ϵ = Inf
sync1_abc = repressilator * aut1_abc
println("After creating sync model")
@show aut1_abc.ϵ
b_sim_aut1_abc = @benchmark (σ = simulate($(sync1_abc)))
@btime (σ = simulate($(sync1_abc)))
println("After bench simulate sync model")
@show aut1_abc.ϵ
@show minimum(b_sim_aut1_abc), mean(b_sim_aut1_abc), maximum(b_sim_aut1_abc)
b_vol_sim_aut1_abc = @benchmark (σ = volatile_simulate($(sync1_abc)))
@btime (σ = volatile_simulate($(sync1_abc)))
println("After bench volatile_simulate sync model")
@show aut1_abc.ϵ
@show minimum(b_vol_sim_aut1_abc), mean(b_vol_sim_aut1_abc), maximum(b_vol_sim_aut1_abc)
#=
println("Memory test")
Profile.clear_malloc_data()
......
......@@ -286,8 +286,10 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym
end
function volatile_simulate(m::$(model_name), A::$(lha_name), p_sim::AbstractVector{Float64},
epsilon::Float64, verbose::Bool)
if $(Meta.quot(lha_name)) == :ABCEuclideanDistanceAutomaton A.ϵ = epsilon end
epsilon::Union{Nothing,Float64}, verbose::Bool)
if $(Meta.quot(lha_name)) == :ABCEuclideanDistanceAutomaton && epsilon != nothing
A.ϵ = epsilon
end
x0 = getfield(m, :x0)
t0 = getfield(m, :t0)
time_bound = getfield(m, :time_bound)
......@@ -341,7 +343,7 @@ function generate_code_synchronized_simulation(model_name::Symbol, lha_name::Sym
if isabsorbing break end
end
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
......@@ -362,7 +364,7 @@ in order to improve performance.
It returns the last state of the simulation `S::StateLHA` not a trajectory `σ::SynchronizedTrajectory`.
"""
function volatile_simulate(product::SynchronizedModel;
p::Union{Nothing,AbstractVector{Float64}} = nothing, epsilon::Float64 = 0.0, verbose::Bool = false)
p::Union{Nothing,AbstractVector{Float64}} = nothing, epsilon::Union{Nothing,Float64} = nothing, verbose::Bool = false)
m = product.m
A = product.automaton
p_sim = getfield(m, :p)
......
......@@ -6,42 +6,64 @@ import Distributions: Uniform
MAKE_SECOND_AUTOMATON_TESTS = false
load_model("SIR")
load_model("repressilator")
load_automaton("euclidean_distance_automaton")
load_automaton("abc_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)
aut1 = create_euclidean_distance_automaton(SIR, tml_obs, y_obs, :I)
y_obs_sir = vectorize(simulate(SIR), :I, tml_obs)
aut1 = create_euclidean_distance_automaton(SIR, tml_obs, y_obs_sir, :I)
sync_SIR = SIR * aut1
set_time_bound!(repressilator, 200.0)
y_obs_repr = vectorize(simulate(repressilator), :P1, tml_obs)
aut1 = create_euclidean_distance_automaton(repressilator, tml_obs, y_obs_repr, :P1)
aut1_abc = create_abc_euclidean_distance_automaton(repressilator, tml_obs, y_obs_repr, :P1)
sync_repressilator = repressilator * aut1
sync_abc_repressilator = repressilator * aut1_abc
test_all = true
nbr_sim = 100
nbr_sim = 10
for i = 1:nbr_sim
let σ, test
σ = simulate(sync_SIR)
test = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
test = euclidean_distance(σ, :I, tml_obs, y_obs_sir) == σ.state_lha_end[:d]
if !test
@show tml_obs
@show euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
@show euclidean_distance(σ, :I, tml_obs, y_obs_sir), σ.state_lha_end[:d]
@show σ
end
σ = simulate(sync_repressilator)
test = test && euclidean_distance(σ, :P1, tml_obs, y_obs_repr) == σ.state_lha_end[:d]
if !test
@show tml_obs
@show euclidean_distance(σ, :P1, tml_obs, y_obs_repr), σ.state_lha_end[:d]
#@show σ
end
σ = simulate(sync_abc_repressilator)
test = test && euclidean_distance(σ, :P1, tml_obs, y_obs_repr) == σ.state_lha_end[:d]
if !test
@show tml_obs
@show euclidean_distance(σ, :P1, tml_obs, y_obs_repr), σ.state_lha_end[:d]
#@show σ
end
global test_all = test_all && test
end
end
if MAKE_SECOND_AUTOMATON_TESTS
aut2 = create_euclidean_distance_automaton_2(SIR, tml_obs, y_obs, :I)
aut2 = create_euclidean_distance_automaton_2(SIR, tml_obs, y_obs_sir, :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]
test2 = euclidean_distance(σ, :I, tml_obs, y_obs_sir) == σ.state_lha_end[:d]
if !test2
@show euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
@show euclidean_distance(σ, :I, tml_obs, y_obs_sir), σ.state_lha_end[:d]
@show σ
end
global test_all = test_all && test2
......@@ -51,35 +73,55 @@ 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)
y_obs_sir = vectorize(simulate(SIR), :I, tml_obs)
aut1 = create_euclidean_distance_automaton(SIR, tml_obs, y_obs_sir, :I)
sync_SIR = SIR * aut1
set_time_bound!(repressilator, 200.0)
y_obs_repr = vectorize(simulate(repressilator), :P1, tml_obs)
aut1 = create_euclidean_distance_automaton(repressilator, tml_obs, y_obs_repr, :P1)
aut1_abc = create_abc_euclidean_distance_automaton(repressilator, tml_obs, y_obs_repr, :P1)
sync_repressilator = repressilator * aut1
sync_abc_repressilator = repressilator * aut1_abc
test_all = true
nbr_sim = 100
nbr_sim = 10
for i = 1:nbr_sim
let σ, test
σ = simulate(sync_SIR)
test = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
test = euclidean_distance(σ, :I, tml_obs, y_obs_sir) == σ.state_lha_end[:d]
if !test
@show tml_obs
@show euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
@show euclidean_distance(σ, :I, tml_obs, y_obs_sir), σ.state_lha_end[:d]
@show σ
end
σ = simulate(sync_repressilator)
test = test && euclidean_distance(σ, :P1, tml_obs, y_obs_repr) == σ.state_lha_end[:d]
if !test
@show tml_obs
@show euclidean_distance(σ, :P1, tml_obs, y_obs_repr), σ.state_lha_end[:d]
#@show σ
end
σ = simulate(sync_abc_repressilator)
test = test && euclidean_distance(σ, :P1, tml_obs, y_obs_repr) == σ.state_lha_end[:d]
if !test
@show tml_obs
@show euclidean_distance(σ, :P1, tml_obs, y_obs_repr), σ.state_lha_end[:d]
#@show σ
end
global test_all = test_all && test
end
end
if MAKE_SECOND_AUTOMATON_TESTS
aut2 = create_euclidean_distance_automaton_2(SIR, tml_obs, y_obs, :I)
aut2 = create_euclidean_distance_automaton_2(SIR, tml_obs, y_obs_sir, :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]
test2 = euclidean_distance(σ, :I, tml_obs, y_obs_sir) == σ.state_lha_end[:d]
if !test2
@show euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
@show euclidean_distance(σ, :I, tml_obs, y_obs_sir), σ.state_lha_end[:d]
@show σ
end
global test_all = test_all && test2
......
Markdown is supported
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