From 9190a75cfaf948d4eaf38162644ce3f31c391fdb Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Thu, 3 Dec 2020 13:05:52 +0100
Subject: [PATCH] add tests and fix warning tests + poisson process
---
algorithms/_utils_abc.jl | 59 ++++++++++------
algorithms/abc_smc.jl | 7 +-
core/common.jl | 5 +-
core/model.jl | 1 +
models/ER.jl | 6 +-
models/SIR_tauleap.jl | 2 +-
models/poisson.jl | 36 ++++++++++
tests/automata/accept_R5.jl | 5 +-
.../automata/read_trajectory_last_state_F.jl | 2 +-
.../automata/read_trajectory_last_state_G.jl | 2 +-
tests/automata/sync_simulate_last_state_F.jl | 2 +-
tests/automata/sync_simulate_last_state_G.jl | 2 +-
tests/automaton_abc/R1.jl | 17 +++--
tests/automaton_abc/distributed_R1.jl | 38 +++++++++++
tests/cosmos/distance_G/ER_R5.jl | 2 +-
tests/dist_lp/dist_case_2.jl | 62 +++++++++--------
tests/dist_lp/dist_case_3.jl | 68 ++++++++++---------
tests/dist_lp/dist_sim_sir.jl | 36 +++++-----
tests/run_abc_smc.jl | 9 +++
tests/run_all.jl | 1 +
tests/run_cosmos.jl | 6 +-
tests/run_unit.jl | 1 +
tests/unit/simulate_available_models.jl | 15 ++++
23 files changed, 260 insertions(+), 124 deletions(-)
create mode 100644 models/poisson.jl
create mode 100644 tests/automaton_abc/distributed_R1.jl
create mode 100644 tests/run_abc_smc.jl
create mode 100644 tests/unit/simulate_available_models.jl
diff --git a/algorithms/_utils_abc.jl b/algorithms/_utils_abc.jl
index ec6e0f6..5c1d373 100644
--- a/algorithms/_utils_abc.jl
+++ b/algorithms/_utils_abc.jl
@@ -1,12 +1,19 @@
function _init_abc_automaton!(mat_p_old::Matrix{Float64}, vec_dist::Vector{Float64},
- pm::ParametricModel, str_var_aut::String)
+ pm::ParametricModel, str_var_aut::String;
+ l_obs::Union{Nothing,AbstractVector} = nothing,
+ func_dist::Union{Nothing,Function} = nothing)
vec_p = zeros(pm.df)
for i = eachindex(vec_dist)
draw!(vec_p, pm)
mat_p_old[:,i] = vec_p
- S = volatile_simulate(pm, vec_p)
- vec_dist[i] = S[str_var_aut]
+ if l_obs == nothing
+ S = volatile_simulate(pm, vec_p)
+ vec_dist[i] = S[str_var_aut]
+ else
+ σ = simulate(pm, vec_p)
+ vec_dist[i] = func_dist(σ, l_obs)
+ end
end
end
@@ -38,23 +45,26 @@ function _task_worker!(d_mat_p::DArray{Float64,2}, d_vec_dist::DArray{Float64,1}
pm::ParametricModel, epsilon::Float64,
wl_old::Vector{Float64}, mat_p_old::Matrix{Float64},
mat_cov::Matrix{Float64}, tree_mat_p::Union{Nothing,KDTree},
- kernel_type::String, str_var_aut::String)
+ kernel_type::String, str_var_aut::String;
+ l_obs::Union{Nothing,AbstractVector} = nothing,
+ func_dist::Union{Nothing,Function} = nothing)
mat_p = localpart(d_mat_p)
vec_dist = localpart(d_vec_dist)
wl_current = localpart(d_wl_current)
l_nbr_sim = zeros(Int, length(vec_dist))
Threads.@threads for i = eachindex(vec_dist)
_update_param!(mat_p, vec_dist, l_nbr_sim, wl_current, i, pm, epsilon,
- wl_old, mat_p_old, mat_cov, tree_mat_p, kernel_type, str_var_aut)
+ wl_old, mat_p_old, mat_cov, tree_mat_p, kernel_type, str_var_aut;
+ l_obs = l_obs, func_dist = func_dist)
end
return sum(l_nbr_sim)
end
function _draw_param_kernel!(vec_p_prime::Vector{Float64},
- vec_p_star::SubArray{Float64,1},
- mat_p_old::Matrix{Float64}, wl_old::Vector{Float64},
- mat_cov::Matrix{Float64},
- tree_mat_p::Union{KDTree,Nothing}, kernel_type::String)
+ vec_p_star::SubArray{Float64,1},
+ mat_p_old::Matrix{Float64}, wl_old::Vector{Float64},
+ mat_cov::Matrix{Float64},
+ tree_mat_p::Union{KDTree,Nothing}, kernel_type::String)
if kernel_type == "mvnormal"
d_mvnormal = MvNormal(vec_p_star, mat_cov)
rand!(d_mvnormal, vec_p_prime)
@@ -71,12 +81,14 @@ function _draw_param_kernel!(vec_p_prime::Vector{Float64},
end
function _update_param!(mat_p::Matrix{Float64}, vec_dist::Vector{Float64},
- l_nbr_sim::Vector{Int}, wl_current::Vector{Float64},
- idx::Int,
- pm::ParametricModel, epsilon::Float64,
- wl_old::Vector{Float64}, mat_p_old::Matrix{Float64},
- mat_cov::Matrix{Float64}, tree_mat_p::Union{Nothing,KDTree},
- kernel_type::String, str_var_aut::String)
+ l_nbr_sim::Vector{Int}, wl_current::Vector{Float64},
+ idx::Int,
+ pm::ParametricModel, epsilon::Float64,
+ wl_old::Vector{Float64}, mat_p_old::Matrix{Float64},
+ mat_cov::Matrix{Float64}, tree_mat_p::Union{Nothing,KDTree},
+ kernel_type::String, str_var_aut::String;
+ l_obs::Union{Nothing,AbstractVector} = nothing,
+ func_dist::Union{Nothing,Function} = nothing)
d_weights = Categorical(wl_old)
dist_sim = Inf
nbr_sim = 0
@@ -88,8 +100,13 @@ function _update_param!(mat_p::Matrix{Float64}, vec_dist::Vector{Float64},
if !insupport(pm, vec_p_prime)
continue
end
- S = volatile_simulate(pm, vec_p_prime)
- dist_sim = S[str_var_aut]
+ if l_obs == nothing
+ S = volatile_simulate(pm, vec_p_prime)
+ dist_sim = S[str_var_aut]
+ else
+ σ = simulate(pm, vec_p)
+ dist_sim = func_dist(σ, l_obs)
+ end
nbr_sim += 1
end
@@ -108,10 +125,10 @@ function _update_param!(mat_p::Matrix{Float64}, vec_dist::Vector{Float64},
end
function _update_weight!(wl_current::Vector{Float64}, idx::Int,
- pm::ParametricModel,
- wl_old::Vector{Float64}, mat_p_old::Matrix{Float64},
- vec_p_prime::Vector{Float64},
- mat_cov_kernel::Matrix{Float64}, kernel_type::String)
+ pm::ParametricModel,
+ wl_old::Vector{Float64}, mat_p_old::Matrix{Float64},
+ vec_p_prime::Vector{Float64},
+ mat_cov_kernel::Matrix{Float64}, kernel_type::String)
denom = 0.0
for j in 1:length(wl_old)
#denom += wl_old[j] * Distributions.pdf(d_normal, inv_sqrt_mat_cov * (vec_p_current - mat_p_old[:,j]))::Float64
diff --git a/algorithms/abc_smc.jl b/algorithms/abc_smc.jl
index adc2426..26cc1d5 100644
--- a/algorithms/abc_smc.jl
+++ b/algorithms/abc_smc.jl
@@ -15,7 +15,7 @@ using Logging
main_pkg_path = get_module_path()
include("$(main_pkg_path)/algorithms/_utils_abc.jl")
-struct ResultAutomatonAbc
+struct ResultAbc
mat_p_end::Matrix{Float64}
mat_cov::Matrix{Float64}
nbr_sim::Int64
@@ -176,7 +176,7 @@ function _abc_smc(pm::ParametricModel, nbr_particles::Int, alpha::Float64,
writedlm(dir_res * "weights_end.csv", wl_old, ',')
writedlm(dir_res * "mat_p_end.csv", mat_p_old, ',')
end
- r = ResultAutomatonAbc(mat_p, mat_cov, nbr_tot_sim, time() - begin_time, vec_dist, epsilon, wl_old, l_ess)
+ r = ResultAbc(mat_p, mat_cov, nbr_tot_sim, time() - begin_time, vec_dist, epsilon, wl_old, l_ess)
return r
end
@@ -277,8 +277,7 @@ function _distributed_abc_smc(pm::ParametricModel, nbr_particles::Int, alpha::Fl
writedlm(dir_res * "weights_end.csv", wl_current, ',')
writedlm(dir_res * "mat_p_end.csv", mat_p, ',')
end
- r = ResultAutomatonAbc(mat_p, mat_cov, nbr_tot_sim, time() - begin_time, convert(Array, d_vec_dist), epsilon, wl_current, l_ess)
+ r = ResultAbc(mat_p, mat_cov, nbr_tot_sim, time() - begin_time, convert(Array, d_vec_dist), epsilon, wl_current, l_ess)
return r
end
-
diff --git a/core/common.jl b/core/common.jl
index 3c87419..1b878b9 100644
--- a/core/common.jl
+++ b/core/common.jl
@@ -100,8 +100,9 @@ function ContinuousTimeModel(d::Int, k::Int, map_var_idx::Dict, map_param_idx::D
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, isabsorbing, time_bound, buffer_size, estim_min_states)
+ new_model = 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, estim_min_states)
+ @assert check_consistency(new_model)
+ return new_model
end
LHA(A::LHA, map_var::Dict{String,Int}) = LHA(A.l_transitions, A.l_loc, A.Λ,
diff --git a/core/model.jl b/core/model.jl
index de7d6ad..9350f3a 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -358,6 +358,7 @@ end
isbounded(m::Model) = get_proba_model(m).time_bound < Inf
function check_consistency(m::ContinuousTimeModel)
@assert m.d == length(m.map_var_idx)
+ @assert m.d == length(m.x0)
@assert m.k == length(m.map_param_idx)
@assert m.k == length(m.p)
@assert length(m.g) <= m.d
diff --git a/models/ER.jl b/models/ER.jl
index da45eba..971c53d 100644
--- a/models/ER.jl
+++ b/models/ER.jl
@@ -3,8 +3,8 @@ import StaticArrays: SVector, SMatrix, @SVector, @SMatrix
d=4
k=3
-dict_var = Dict("E" => 1, "S" => 2, "ES" => 3, "P" => 4)
-dict_p = Dict("k1" => 1, "k2" => 2, "k3" => 3)
+dict_var_ER = Dict("E" => 1, "S" => 2, "ES" => 3, "P" => 4)
+dict_p_ER = Dict("k1" => 1, "k2" => 2, "k3" => 3)
l_tr_ER = ["R1","R2","R3"]
p_ER = [1.0, 1.0, 1.0]
x0_ER = [100, 100, 0, 0]
@@ -48,7 +48,7 @@ isabsorbing_ER(p::Vector{Float64},xn::Vector{Int}) =
(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!,isabsorbing_ER; g=g_ER)
+ER = ContinuousTimeModel(d,k,dict_var_ER,dict_p_ER,l_tr_ER,p_ER,x0_ER,t0_ER,ER_f!,isabsorbing_ER; g=g_ER)
function create_ER(new_p::Vector{Float64})
ER_new = deepcopy(ER)
diff --git a/models/SIR_tauleap.jl b/models/SIR_tauleap.jl
index 1ed78f2..b5dfd51 100644
--- a/models/SIR_tauleap.jl
+++ b/models/SIR_tauleap.jl
@@ -3,7 +3,7 @@ import StaticArrays: SVector, SMatrix, @SVector, @SMatrix
import Distributions: Poisson, rand
d=3
-k=2
+k=3
dict_var_SIR_tauleap = Dict("S" => 1, "I" => 2, "R" => 3)
dict_p_SIR_tauleap = Dict("ki" => 1, "kr" => 2, "tau" => 3)
l_tr_SIR_tauleap = ["U"]
diff --git a/models/poisson.jl b/models/poisson.jl
new file mode 100644
index 0000000..f5991ae
--- /dev/null
+++ b/models/poisson.jl
@@ -0,0 +1,36 @@
+
+import StaticArrays: SVector, SMatrix, @SVector, @SMatrix
+import Distributions: Poisson, rand
+
+d=1
+k=1
+dict_var_poisson = Dict("N" => 1)
+dict_p_poisson = Dict("λ" => 1)
+l_tr_poisson = ["R"]
+p_poisson = [5.0]
+x0_poisson = [0]
+t0_poisson = 0.0
+function poisson_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}},
+ xn::Vector{Int}, tn::Float64, p::Vector{Float64})
+
+ u1 = rand()
+ tau = (-log(u1)/p[1])
+ xnplus1[1] += 1
+ l_t[1] = tn + tau
+ l_tr[1] = "U"
+end
+isabsorbing_poisson(p::Vector{Float64}, xn::Vector{Int}) = p[1] === 0.0
+g_poisson = ["N"]
+
+poisson = ContinuousTimeModel(d,k,dict_var_poisson,dict_p_poisson,l_tr_poisson,
+ p_poisson,x0_poisson,t0_poisson,
+ poisson_f!,isabsorbing_poisson; g=g_poisson, time_bound=1.0)
+function create_poisson(new_p::Vector{Float64})
+ poisson_new = deepcopy(poisson)
+ @assert length(poisson_new.p) == length(new_p)
+ set_param!(poisson_new, new_p)
+ return poisson_new
+end
+
+export poisson, create_poisson
+
diff --git a/tests/automata/accept_R5.jl b/tests/automata/accept_R5.jl
index a9ab3fb..08edfa1 100644
--- a/tests/automata/accept_R5.jl
+++ b/tests/automata/accept_R5.jl
@@ -16,11 +16,12 @@ set_param!(ER, "k2", 40.0)
test_all = true
sync_ER = ER*A_G
-σ = nothing
nb_sim = 1000
for i = 1:nb_sim
- global σ = simulate(sync_ER)
+ local σ = simulate(sync_ER)
+ global test_all = test_all && isaccepted(σ)
if !isaccepted(σ)
+ @show σ
error("Ouch")
end
end
diff --git a/tests/automata/read_trajectory_last_state_F.jl b/tests/automata/read_trajectory_last_state_F.jl
index c1c799f..2ad39be 100644
--- a/tests/automata/read_trajectory_last_state_F.jl
+++ b/tests/automata/read_trajectory_last_state_F.jl
@@ -24,7 +24,7 @@ end
test_all = true
nbr_sim = 10000
for i = 1:nbr_sim
- test = test_last_state(A_F, SIR)
+ local test = test_last_state(A_F, SIR)
global test_all = test_all && test
end
diff --git a/tests/automata/read_trajectory_last_state_G.jl b/tests/automata/read_trajectory_last_state_G.jl
index cdd5a43..8578c0a 100644
--- a/tests/automata/read_trajectory_last_state_G.jl
+++ b/tests/automata/read_trajectory_last_state_G.jl
@@ -24,7 +24,7 @@ end
test_all = true
nbr_sim = 10000
for i = 1:nbr_sim
- test = test_last_state(A_G, ER)
+ local test = test_last_state(A_G, ER)
global test_all = test_all && test
end
diff --git a/tests/automata/sync_simulate_last_state_F.jl b/tests/automata/sync_simulate_last_state_F.jl
index 5a1cf7a..0b4ae2b 100644
--- a/tests/automata/sync_simulate_last_state_F.jl
+++ b/tests/automata/sync_simulate_last_state_F.jl
@@ -28,7 +28,7 @@ end
test_all = true
nbr_sim = 10000
for i = 1:nbr_sim
- test = test_last_state(A_F, sync_SIR)
+ local test = test_last_state(A_F, sync_SIR)
global test_all = test_all && test
end
diff --git a/tests/automata/sync_simulate_last_state_G.jl b/tests/automata/sync_simulate_last_state_G.jl
index 1a95edf..2a891db 100644
--- a/tests/automata/sync_simulate_last_state_G.jl
+++ b/tests/automata/sync_simulate_last_state_G.jl
@@ -18,7 +18,7 @@ end
test_all = true
nbr_sim = 10000
for i = 1:nbr_sim
- test = test_last_state(A_G, sync_ER)
+ local test = test_last_state(A_G, sync_ER)
global test_all = test_all && test
end
diff --git a/tests/automaton_abc/R1.jl b/tests/automaton_abc/R1.jl
index 6bf0707..f980407 100644
--- a/tests/automaton_abc/R1.jl
+++ b/tests/automaton_abc/R1.jl
@@ -1,15 +1,24 @@
using MarkovProcesses
-using Plots
load_model("ER")
load_automaton("automaton_F")
A_F_R1 = create_automaton_F(ER, 50.0, 75.0, 0.025, 0.05, "P")
sync_ER = A_F_R1*ER
pm_sync_ER = ParametricModel(sync_ER, ("k3", Uniform(0.0, 100.0)))
+nbr_pa = 502
-r = automaton_abc(pm_sync_ER; nbr_particles = 1000)
+r = automaton_abc(pm_sync_ER; nbr_particles = nbr_pa)
-histogram(r.mat_p_end', weights = r.weights, normalize = :density)
-png("R1_hist.png")
+test = size(r.mat_p_end)[1] == pm_sync_ER.df &&
+ size(r.mat_p_end)[2] == nbr_pa &&
+ length(r.vec_dist) == nbr_pa &&
+ length(r.weights) == nbr_pa &&
+ r.epsilon == 0.0
+
+return test
+
+#using Plots
+#histogram(r.mat_p_end', weights = r.weights, normalize = :density)
+#png("R1_hist.png")
diff --git a/tests/automaton_abc/distributed_R1.jl b/tests/automaton_abc/distributed_R1.jl
new file mode 100644
index 0000000..583a80c
--- /dev/null
+++ b/tests/automaton_abc/distributed_R1.jl
@@ -0,0 +1,38 @@
+
+using Distributed
+using MarkovProcesses
+begin_procs = nprocs()
+if begin_procs == 1
+ addprocs(2)
+end
+path_module = get_module_path() * "/core"
+
+@everywhere begin
+ path_module = $(path_module)
+ push!(LOAD_PATH, path_module)
+ using MarkovProcesses
+ load_model("ER")
+ load_automaton("automaton_F")
+end
+A_F_R1 = create_automaton_F(ER, 50.0, 75.0, 0.025, 0.05, "P")
+sync_ER = A_F_R1*ER
+pm_sync_ER = ParametricModel(sync_ER, ("k3", Uniform(0.0, 100.0)))
+nbr_pa = 404
+
+r = automaton_abc(pm_sync_ER; nbr_particles = nbr_pa)
+
+if begin_procs == 1
+ rmprocs(workers())
+end
+
+test = size(r.mat_p_end)[1] == pm_sync_ER.df &&
+ size(r.mat_p_end)[2] == nbr_pa &&
+ length(r.vec_dist) == nbr_pa &&
+ length(r.weights) == nbr_pa &&
+ r.epsilon == 0.0
+
+return test
+
+#histogram(r.mat_p_end', weights = r.weights, normalize = :density)
+#png("R1_hist.png")
+
diff --git a/tests/cosmos/distance_G/ER_R5.jl b/tests/cosmos/distance_G/ER_R5.jl
index c62185c..dfcf423 100644
--- a/tests/cosmos/distance_G/ER_R5.jl
+++ b/tests/cosmos/distance_G/ER_R5.jl
@@ -10,7 +10,7 @@
observe_all!(ER)
ER.buffer_size = 100
load_automaton("automaton_G")
- width = 0.5
+ width = 0.2
level = 0.95
x1, x2, t1, t2 = 50.0, 100.0, 0.0, 0.8
A_G = create_automaton_G(model, x1, x2, t1, t2, "E")
diff --git a/tests/dist_lp/dist_case_2.jl b/tests/dist_lp/dist_case_2.jl
index 102280b..1bdff98 100644
--- a/tests/dist_lp/dist_case_2.jl
+++ b/tests/dist_lp/dist_case_2.jl
@@ -6,36 +6,38 @@ load_model("SIR")
test_all = true
p=2
for p = 1:2
- res = (4-1)^p * 0.5 + (4-2)^p * 0.1 + 0 + (5-3)^p * 0.1 + (5-1)^p * (1.4 - 0.8) + (3-1)^p * (3.0-1.4)
- res = res^(1/p)
-
- y_obs = [1, 2, 5, 3, 3]
- t_y = [0.0, 0.5, 0.7, 1.4, 3.0]
- values = [zeros(length(y_obs))]
- values[1] = y_obs
- l_tr = Vector{Nothing}(nothing, length(y_obs))
- σ2 = Trajectory(SIR, values, t_y, l_tr)
-
- x_obs = [4, 2, 3, 1, 1]
- t_x = [0.0, 0.6, 0.7, 0.8, 3.0]
- values = [zeros(length(x_obs))]
- values[1] = x_obs
- l_tr = Vector{Nothing}(nothing, length(x_obs))
- σ1 = Trajectory(SIR, values, t_x, l_tr)
-
- test_1 = isapprox(dist_lp(x_obs, t_x, y_obs, t_y; p=p), res; atol = 1E-10)
- test_1 = test_1 && isapprox(dist_lp(σ1,σ2; p=p), res; atol = 1E-10)
- test_1_bis = isapprox(dist_lp(y_obs, t_y, x_obs, t_x; p=p), res; atol = 1E-10)
- test_1_bis = test_1_bis && isapprox(dist_lp(σ2,σ1; p=p), res; atol = 1E-10)
-
- f_x(t::Real) = MarkovProcesses._f_step(x_obs, t_x, t)
- f_y(t::Real) = MarkovProcesses._f_step(y_obs, t_y, t)
- diff_f(t) = abs(f_x(t) - f_y(t))^p
- int, err = quadgk(diff_f, 0.0, 3.0)
- res_int = int^(1/p)
-
- test_2 = isapprox(res, res_int; atol = err)
- global test_all = test_all && test_1 && test_1_bis && test_2
+ let x_obs, y_obs, t_x, t_y, σ1, σ2, test_1, test_1_bis, test_2
+ res = (4-1)^p * 0.5 + (4-2)^p * 0.1 + 0 + (5-3)^p * 0.1 + (5-1)^p * (1.4 - 0.8) + (3-1)^p * (3.0-1.4)
+ res = res^(1/p)
+
+ y_obs = [1, 2, 5, 3, 3]
+ t_y = [0.0, 0.5, 0.7, 1.4, 3.0]
+ values = [zeros(length(y_obs))]
+ values[1] = y_obs
+ l_tr = Vector{Nothing}(nothing, length(y_obs))
+ σ2 = Trajectory(SIR, values, t_y, l_tr)
+
+ x_obs = [4, 2, 3, 1, 1]
+ t_x = [0.0, 0.6, 0.7, 0.8, 3.0]
+ values = [zeros(length(x_obs))]
+ values[1] = x_obs
+ l_tr = Vector{Nothing}(nothing, length(x_obs))
+ σ1 = Trajectory(SIR, values, t_x, l_tr)
+
+ test_1 = isapprox(dist_lp(x_obs, t_x, y_obs, t_y; p=p), res; atol = 1E-10)
+ test_1 = test_1 && isapprox(dist_lp(σ1,σ2; p=p), res; atol = 1E-10)
+ test_1_bis = isapprox(dist_lp(y_obs, t_y, x_obs, t_x; p=p), res; atol = 1E-10)
+ test_1_bis = test_1_bis && isapprox(dist_lp(σ2,σ1; p=p), res; atol = 1E-10)
+
+ f_x(t::Real) = MarkovProcesses._f_step(x_obs, t_x, t)
+ f_y(t::Real) = MarkovProcesses._f_step(y_obs, t_y, t)
+ diff_f(t) = abs(f_x(t) - f_y(t))^p
+ int, err = quadgk(diff_f, 0.0, 3.0)
+ res_int = int^(1/p)
+
+ test_2 = isapprox(res, res_int; atol = err)
+ global test_all = test_all && test_1 && test_1_bis && test_2
+ end
end
return test_all
diff --git a/tests/dist_lp/dist_case_3.jl b/tests/dist_lp/dist_case_3.jl
index 8254be0..f22dabf 100644
--- a/tests/dist_lp/dist_case_3.jl
+++ b/tests/dist_lp/dist_case_3.jl
@@ -5,39 +5,41 @@ load_model("SIR")
test_all = true
for p = 1:2
- x_obs =[5, 6, 5, 4, 3, 2, 1, 1]
- t_x = [0.0, 3.10807, 4.29827, 4.40704, 5.67024, 7.1299, 11.2763, 20.0]
- values = [zeros(length(x_obs))]
- values[1] = x_obs
- l_tr = Vector{Nothing}(nothing, length(x_obs))
- σ1 = Trajectory(SIR, values, t_x, l_tr)
-
- y_obs =[5, 4, 5, 4, 3, 4, 3, 2, 3, 2, 1, 2, 3, 4, 3, 4, 4]
- t_y = [0.0, 0.334082, 1.21012, 1.40991, 1.58866, 2.45879, 2.94545, 4.66746, 5.44723, 5.88066, 7.25626, 11.4036, 13.8373, 17.1363, 17.8193, 18.7613, 20.0]
- values = [zeros(length(y_obs))]
- values[1] = y_obs
- l_tr = Vector{Nothing}(nothing, length(y_obs))
- σ2 = Trajectory(SIR, values, t_y, l_tr)
-
- f_x(t::Real) = MarkovProcesses._f_step(x_obs, t_x, t)
- f_y(t::Real) = MarkovProcesses._f_step(y_obs, t_y, t)
- diff_f(t) = abs(f_x(t) - f_y(t))^p
- int, err = quadgk(diff_f, 0.0, 20.0, rtol=1e-10)
- int_riemann = MarkovProcesses._riemann_sum(diff_f, 0.0, 20.0, 1E-5)
- int_riemann = int_riemann^(1/p)
-
- res1 = dist_lp(x_obs, t_x, y_obs, t_y; p=p)
- res2 = dist_lp(σ1,σ2; p=p)
- res1_bis = dist_lp(y_obs, t_y, x_obs, t_x; p=p)
- res2_bis = dist_lp(σ2,σ1; p=p)
- test_1 = isapprox(res1, int_riemann; atol = 1E-3)
- test_1 = test_1 && isapprox(res2, int_riemann; atol = 1E-3)
- test_1_bis = isapprox(res1_bis, int_riemann; atol = 1E-3)
- test_1_bis = test_1_bis && isapprox(res2_bis, int_riemann; atol = 1E-3)
-
- test_2 = res1 == res2 == res1_bis == res2_bis
-
- global test_all = test_all && test_1 && test_1_bis && test_2
+ let x_obs, y_obs, t_x, t_y, σ1, σ2, test_1, test_1_bis, test_2
+ x_obs =[5, 6, 5, 4, 3, 2, 1, 1]
+ t_x = [0.0, 3.10807, 4.29827, 4.40704, 5.67024, 7.1299, 11.2763, 20.0]
+ values = [zeros(length(x_obs))]
+ values[1] = x_obs
+ l_tr = Vector{Nothing}(nothing, length(x_obs))
+ σ1 = Trajectory(SIR, values, t_x, l_tr)
+
+ y_obs =[5, 4, 5, 4, 3, 4, 3, 2, 3, 2, 1, 2, 3, 4, 3, 4, 4]
+ t_y = [0.0, 0.334082, 1.21012, 1.40991, 1.58866, 2.45879, 2.94545, 4.66746, 5.44723, 5.88066, 7.25626, 11.4036, 13.8373, 17.1363, 17.8193, 18.7613, 20.0]
+ values = [zeros(length(y_obs))]
+ values[1] = y_obs
+ l_tr = Vector{Nothing}(nothing, length(y_obs))
+ σ2 = Trajectory(SIR, values, t_y, l_tr)
+
+ f_x(t::Real) = MarkovProcesses._f_step(x_obs, t_x, t)
+ f_y(t::Real) = MarkovProcesses._f_step(y_obs, t_y, t)
+ diff_f(t) = abs(f_x(t) - f_y(t))^p
+ int, err = quadgk(diff_f, 0.0, 20.0, rtol=1e-10)
+ int_riemann = MarkovProcesses._riemann_sum(diff_f, 0.0, 20.0, 1E-5)
+ int_riemann = int_riemann^(1/p)
+
+ res1 = dist_lp(x_obs, t_x, y_obs, t_y; p=p)
+ res2 = dist_lp(σ1,σ2; p=p)
+ res1_bis = dist_lp(y_obs, t_y, x_obs, t_x; p=p)
+ res2_bis = dist_lp(σ2,σ1; p=p)
+ test_1 = isapprox(res1, int_riemann; atol = 1E-3)
+ test_1 = test_1 && isapprox(res2, int_riemann; atol = 1E-3)
+ test_1_bis = isapprox(res1_bis, int_riemann; atol = 1E-3)
+ test_1_bis = test_1_bis && isapprox(res2_bis, int_riemann; atol = 1E-3)
+
+ test_2 = res1 == res2 == res1_bis == res2_bis
+
+ global test_all = test_all && test_1 && test_1_bis && test_2
+ end
end
return test_all
diff --git a/tests/dist_lp/dist_sim_sir.jl b/tests/dist_lp/dist_sim_sir.jl
index 13e3889..61008a3 100644
--- a/tests/dist_lp/dist_sim_sir.jl
+++ b/tests/dist_lp/dist_sim_sir.jl
@@ -9,23 +9,25 @@ test_all = true
for p = 1:2
nb_sim = 10
for i = 1:nb_sim
- σ1 = simulate(SIR)
- σ2 = simulate(SIR)
- d = dist_lp(σ1, σ2, "I"; p = p)
- d2 = dist_lp(σ1, σ2; p = p)
-
- f_x(t::Float64) = MarkovProcesses._f_step(σ1["I"], times(σ1), t)
- f_y(t::Float64) = MarkovProcesses._f_step(σ2["I"], times(σ2), t)
- diff_f(t) = abs(f_x(t) - f_y(t))^p
- int_riemann = MarkovProcesses._riemann_sum(diff_f, SIR.t0, SIR.time_bound, 1E-3)
- res_int_riemann = int_riemann^(1/p)
-
- test = isapprox(d, res_int_riemann; atol = 1E-1)
- test2 = isapprox(d2, res_int_riemann; atol = 1E-1)
- #@show d, res_int_riemann
- #@show test
-
- global test_all = test_all && test && test2
+ let σ1, σ2, test, test2
+ σ1 = simulate(SIR)
+ σ2 = simulate(SIR)
+ d = dist_lp(σ1, σ2, "I"; p = p)
+ d2 = dist_lp(σ1, σ2; p = p)
+
+ f_x(t::Float64) = MarkovProcesses._f_step(σ1["I"], times(σ1), t)
+ f_y(t::Float64) = MarkovProcesses._f_step(σ2["I"], times(σ2), t)
+ diff_f(t) = abs(f_x(t) - f_y(t))^p
+ int_riemann = MarkovProcesses._riemann_sum(diff_f, SIR.t0, SIR.time_bound, 1E-3)
+ res_int_riemann = int_riemann^(1/p)
+
+ test = isapprox(d, res_int_riemann; atol = 1E-1)
+ test2 = isapprox(d2, res_int_riemann; atol = 1E-1)
+ #@show d, res_int_riemann
+ #@show test
+
+ global test_all = test_all && test && test2
+ end
end
end
diff --git a/tests/run_abc_smc.jl b/tests/run_abc_smc.jl
new file mode 100644
index 0000000..c922840
--- /dev/null
+++ b/tests/run_abc_smc.jl
@@ -0,0 +1,9 @@
+
+using Test
+
+@testset "ABC SMC and automaton-ABC tests" begin
+ @test include("automaton_abc/R1.jl")
+ #test_distributed_R1 = include("automaton_abc/distributed_R1.jl")
+ #@test test_distributed_R1
+end
+
diff --git a/tests/run_all.jl b/tests/run_all.jl
index 67ed817..544afd7 100644
--- a/tests/run_all.jl
+++ b/tests/run_all.jl
@@ -4,4 +4,5 @@ include("run_simulation.jl")
include("run_dist_lp.jl")
include("run_automata.jl")
include("run_cosmos.jl")
+include("run_abc_smc.jl")
diff --git a/tests/run_cosmos.jl b/tests/run_cosmos.jl
index bc2ad79..b42c650 100644
--- a/tests/run_cosmos.jl
+++ b/tests/run_cosmos.jl
@@ -2,7 +2,9 @@
using Test
@testset "Cosmos tests" begin
- @test include("cosmos/distance_F/ER_1D.jl")
- @test include("cosmos/distance_G/ER_R5.jl")
+ test_ER_1D = include("cosmos/distance_F/ER_1D.jl")
+ test_ER_R5 = include("cosmos/distance_G/ER_R5.jl")
+ @test test_ER_1D
+ @test test_ER_R5
end
diff --git a/tests/run_unit.jl b/tests/run_unit.jl
index 1fbcbf6..24cd7f7 100644
--- a/tests/run_unit.jl
+++ b/tests/run_unit.jl
@@ -33,6 +33,7 @@ using Test
@test include("unit/set_param.jl")
@test include("unit/side_effects_1.jl")
+ @test include("unit/simulate_available_models.jl")
@test include("unit/simulate_sir.jl")
@test include("unit/simulate_sir_bounded.jl")
@test include("unit/simulate_er.jl")
diff --git a/tests/unit/simulate_available_models.jl b/tests/unit/simulate_available_models.jl
new file mode 100644
index 0000000..8c6bbec
--- /dev/null
+++ b/tests/unit/simulate_available_models.jl
@@ -0,0 +1,15 @@
+
+using MarkovProcesses
+
+load_model("SIR")
+load_model("SIR_tauleap")
+load_model("ER")
+load_model("poisson")
+
+simulate(SIR)
+simulate(ER)
+simulate(SIR_tauleap)
+simulate(poisson)
+
+return true
+
--
GitLab