From 119c1c41de69deb63c5091e196f0e71b6bc4cb39 Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Tue, 16 Feb 2021 12:02:40 +0100
Subject: [PATCH] add of a new automaton for computing euclidean distances +
tests
---
automata/euclidean_distance_automaton_2.jl | 119 ++++++++++++++++++++
tests/automata/euclidean_distance.jl | 52 ++++++++-
tests/automata/euclidean_distance_single.jl | 11 +-
3 files changed, 175 insertions(+), 7 deletions(-)
create mode 100644 automata/euclidean_distance_automaton_2.jl
diff --git a/automata/euclidean_distance_automaton_2.jl b/automata/euclidean_distance_automaton_2.jl
new file mode 100644
index 0000000..ed69d69
--- /dev/null
+++ b/automata/euclidean_distance_automaton_2.jl
@@ -0,0 +1,119 @@
+
+function create_euclidean_distance_automaton_2(m::ContinuousTimeModel, timeline::AbstractVector{Float64}, observations::AbstractVector{Float64}, sym_obs::VariableModel)
+ # Requirements for the automaton
+ @assert sym_obs in m.g "$(sym_obs) is not observed."
+ @assert length(timeline) == length(observations) "Timeline and observations vectors don't have the same length"
+ nbr_observations = length(observations)
+
+ # Locations
+ locations = [:l0, :lfinal]
+ for i = 1:nbr_observations
+ push!(locations, Symbol("l$(i)"))
+ end
+
+ ## Invariant predicates
+ @everywhere true_inv_predicate(x::Vector{Int}) = true
+ Λ_F = Dict{Location, Function}()
+ for loc in locations
+ Λ_F[loc] = getfield(Main, :true_inv_predicate)
+ end
+
+ ## Init and final loc
+ locations_init = [:l0]
+ locations_final = [:lfinal]
+
+ map_var_automaton_idx = Dict{VariableAutomaton,Int}(:t => 1, :n => 2, :d => 3)
+ vector_flow = [1.0, 0.0, 0.0]
+ flow = Dict{Location, Vector{Float64}}()
+ for loc in locations
+ flow[loc] = vector_flow
+ end
+
+ ## Edges
+ map_edges = Dict{Location, Dict{Location, Vector{Edge}}}()
+ for loc in locations
+ map_edges[loc] = Dict{Location, Vector{Edge}}()
+ end
+
+ idx_obs_var = getfield(m, :map_var_idx)[sym_obs]
+ to_idx(var::Symbol) = map_var_automaton_idx[var]
+ nbr_rand = rand(1:1000)
+ basename_func = "$(replace(m.name, ' '=>'_'))_$(nbr_rand)"
+ basename_func = replace(basename_func, '-'=>'_')
+ func_name(type_func::Symbol, from_loc::Location, to_loc::Location, edge_number::Int) =
+ Symbol("$(type_func)_eucl_dist_aut_2_$(basename_func)_$(from_loc)$(to_loc)_$(edge_number)$(type_func == :us ? "!" : "")")
+ loc_nbr_obs = Symbol("l$(nbr_observations)")
+ meta_elementary_functions = quote
+ # l0 loc
+ # l0 => l1
+ @everywhere $(func_name(:cc, :l0, :l1, 1))(S::StateLHA, x::Vector{Int}, p::Vector{Float64}) = true
+ @everywhere $(func_name(:us, :l0, :l1, 1))(S::StateLHA, x::Vector{Int}, p::Vector{Float64}) =
+ (setfield!(S, :loc, Symbol("l1"));
+ setindex!(getfield(S, :values), x[$(idx_obs_var)], $(to_idx(:n)));
+ setindex!(getfield(S, :values), 0.0, $(to_idx(:d))))
+
+ # lnbr_obs => lfinal
+ @everywhere $(func_name(:cc, loc_nbr_obs, :lfinal, 1))(S::StateLHA, x::Vector{Int}, p::Vector{Float64}) =
+ get_value(S, $(to_idx(:t))) >= $(timeline[nbr_observations])
+ @everywhere $(func_name(:us, loc_nbr_obs, :lfinal, 1))(S::StateLHA, x::Vector{Int}, p::Vector{Float64}) =
+ (setfield!(S, :loc, Symbol("lfinal"));
+ setindex!(getfield(S, :values), get_value(S, $(to_idx(:d))) + (get_value(S, $(to_idx(:n)))-$(observations[nbr_observations]))^2,
+ $(to_idx(:d)));
+ setindex!(getfield(S, :values), sqrt(get_value(S, $(to_idx(:d)))), $(to_idx(:d))))
+
+ # lnbr_obs => lnbr_obs
+ @everywhere $(func_name(:cc, loc_nbr_obs, loc_nbr_obs, 1))(S::StateLHA, x::Vector{Int}, p::Vector{Float64}) = true
+ @everywhere $(func_name(:us, loc_nbr_obs, loc_nbr_obs, 1))(S::StateLHA, x::Vector{Int}, p::Vector{Float64}) =
+ (setindex!(getfield(S, :values), x[$(idx_obs_var)], $(to_idx(:n))))
+ end
+ eval(meta_elementary_functions)
+ # l0 loc
+ # l0 => l1
+ edge1 = Edge([nothing], getfield(Main, func_name(:cc, :l0, :l1, 1)), getfield(Main, func_name(:us, :l0, :l1, 1)))
+ map_edges[:l0][:l1] = [edge1]
+
+ # lnbr_obs => lfinal
+ edge1 = Edge([nothing], getfield(Main, func_name(:cc, loc_nbr_obs, :lfinal, 1)), getfield(Main, func_name(:us, loc_nbr_obs, :lfinal, 1)))
+ map_edges[loc_nbr_obs][:lfinal] = [edge1]
+ # lnbr_obs => lnbr_obs
+ edge1 = Edge([:ALL], getfield(Main, func_name(:cc, loc_nbr_obs, loc_nbr_obs, 1)), getfield(Main, func_name(:us, loc_nbr_obs, loc_nbr_obs, 1)))
+ map_edges[loc_nbr_obs][loc_nbr_obs] = [edge1]
+
+ for i = 1:(nbr_observations-1)
+ loci = Symbol("l$(i)")
+ locip1 = Symbol("l$(i+1)")
+ meta_elementary_functions_loci = quote
+ # l1 loc
+ # l1 => l1
+ # Defined below
+ @everywhere $(func_name(:cc, loci, locip1, 1))(S::StateLHA, x::Vector{Int}, p::Vector{Float64}) =
+ get_value(S, $(to_idx(:t))) >= $(timeline[i])
+ @everywhere $(func_name(:us, loci, locip1, 1))(S::StateLHA, x::Vector{Int}, p::Vector{Float64}) =
+ (setfield!(S, :loc, $(Meta.quot(locip1)));
+ setindex!(getfield(S, :values), get_value(S, $(to_idx(:d))) + (get_value(S, $(to_idx(:n)))-$(observations[i]))^2,
+ $(to_idx(:d))))
+
+ @everywhere $(func_name(:cc, loci, loci, 1))(S::StateLHA, x::Vector{Int}, p::Vector{Float64}) = true
+ @everywhere $(func_name(:us, loci, loci, 1))(S::StateLHA, x::Vector{Int}, p::Vector{Float64}) =
+ (setindex!(getfield(S, :values), x[$(idx_obs_var)], $(to_idx(:n))))
+ end
+ eval(meta_elementary_functions_loci)
+
+ # loci => loci+1
+ edge1 = Edge([nothing], getfield(Main, func_name(:cc, loci, locip1, 1)), getfield(Main, func_name(:us, loci, locip1, 1)))
+ map_edges[loci][locip1] = [edge1]
+ # loci => loci
+ edge1 = Edge([:ALL], getfield(Main, func_name(:cc, loci, loci, 1)), getfield(Main, func_name(:us, loci, loci, 1)))
+ map_edges[loci][loci] = [edge1]
+ end
+
+ ## Constants
+ constants = Dict{Symbol,Float64}(:nbr_obs => nbr_observations)
+
+ A = LHA("Euclidean distance", m.transitions, locations, Λ_F, locations_init, locations_final,
+ map_var_automaton_idx, flow, map_edges, constants, m.map_var_idx)
+ return A
+end
+
+export create_euclidean_distance_automaton_2
+
diff --git a/tests/automata/euclidean_distance.jl b/tests/automata/euclidean_distance.jl
index 2f0cc1f..bc1eff5 100644
--- a/tests/automata/euclidean_distance.jl
+++ b/tests/automata/euclidean_distance.jl
@@ -4,37 +4,77 @@ import LinearAlgebra: dot
import Distributions: Uniform
load_automaton("euclidean_distance_automaton")
+load_automaton("euclidean_distance_automaton_2")
load_model("SIR")
load_model("ER")
+observe_all!(SIR)
+observe_all!(ER)
test_all = true
# SIR model
-nbr_sim = 10
+nbr_sim = 20
for i = 1:nbr_sim
set_param!(SIR, [:ki, :kr], [rand(Uniform(5E-5, 3E-3)), rand(Uniform(5E-3, 0.2))])
- let tml_obs, y_obs, sync_SIR, σ, test
+ let tml_obs, y_obs, sync_SIR, σ, test, test2
tml_obs = rand(Uniform(0.0, 5.0)):1.0:rand(Uniform(50.0, 100.0))
y_obs = vectorize(simulate(SIR), :I, tml_obs)
sync_SIR = SIR * create_euclidean_distance_automaton(SIR, tml_obs, y_obs, :I)
σ = simulate(sync_SIR)
test = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
- #@show test, euclidean_distance(σ, tml_obs, y_obs, :I), σ.state_lha_end[:d]
- global test_all = test_all && test
+ if !test
+ @show test, 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
+ 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
+ global test_all = test_all && test && test2
end
end
# ER model
for i = 1:nbr_sim
- let tml_obs, y_obs, sync_SIR, σ, test
+ let tml_obs, y_obs, sync_SIR, σ, test, test2
set_param!(ER, :k3, rand(Uniform(0.0, 100.0)))
tml_obs = rand(Uniform(0.0, 0.2)):1.0:rand(Uniform(0.5,10.0))
y_obs = vectorize(simulate(ER), :P, tml_obs)
sync_ER = ER * create_euclidean_distance_automaton(ER, tml_obs, y_obs, :P)
σ = simulate(sync_ER)
test = euclidean_distance(σ, :P, tml_obs, y_obs) == σ.state_lha_end[:d]
+ if !test
+ @show test, 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
+ 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
#@show test, euclidean_distance(σ, tml_obs, y_obs, :P), σ.state_lha_end[:d]
- global test_all = test_all && test
+ global test_all = test_all && test && test2
end
end
diff --git a/tests/automata/euclidean_distance_single.jl b/tests/automata/euclidean_distance_single.jl
index f1923c5..fdebb5a 100644
--- a/tests/automata/euclidean_distance_single.jl
+++ b/tests/automata/euclidean_distance_single.jl
@@ -4,12 +4,21 @@ import LinearAlgebra: dot
import Distributions: Uniform
load_automaton("euclidean_distance_automaton")
+load_automaton("euclidean_distance_automaton_2")
load_model("SIR")
tml_obs = 0:0.5:200
set_time_bound!(SIR, 200.0)
y_obs = vectorize(simulate(SIR), :I, tml_obs)
-sync_SIR = SIR * create_euclidean_distance_automaton(SIR, tml_obs, y_obs, :I)
+
+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)
+@show euclidean_distance(σ, :I, tml_obs, y_obs), σ.state_lha_end[:d]
test = euclidean_distance(σ, :I, tml_obs, y_obs) == σ.state_lha_end[:d]
return test
--
GitLab