add tests and fix warning tests + poisson process

algorithms/_utils_abc.jl

 
 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 

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
 
-

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.Λ, 

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

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)

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"]

models/poisson.jl

+
+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
+

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

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
 

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
 

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
 

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
 

tests/automaton_abc/R1.jl

 
 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")
 

tests/automaton_abc/distributed_R1.jl

+
+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")
+

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")  

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

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

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
 

tests/run_abc_smc.jl

+
+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
+

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")