diff --git a/automata/euclidean_distance_automaton.jl b/automata/euclidean_distance_automaton.jl
new file mode 100644
index 0000000000000000000000000000000000000000..703b2224da28a4036728c815bca322140939ad88
--- /dev/null
+++ b/automata/euclidean_distance_automaton.jl
@@ -0,0 +1,108 @@
+
+# l0 loc
+# l0 => l1
+@everywhere cc_eucl_dist_aut_l0l1_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) = true
+
+# l1 loc
+# l1 => l1
+@everywhere cc_eucl_dist_aut_l1l1_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:t] >= constants[Symbol("tml_$(S[:idx])")]
+@everywhere us_eucl_dist_aut_l1l1_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S[:d] += (S[:n] - constants[Symbol("y_$(S[:idx])")])^2;
+ S[:idx] += 1.0)
+
+# l1 => l2
+@everywhere cc_eucl_dist_aut_l1l2_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:idx] >= constants[:nbr_obs] + 1
+@everywhere us_eucl_dist_aut_l1l2_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S.loc = :l2;
+ S[:d] = sqrt(S[:d]))
+
+@everywhere cc_eucl_dist_aut_l1l1_2(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) = true
+
+function create_euclidean_distance_automaton(m::ContinuousTimeModel, timeline::Vector{Float64}, observations::Vector{Float64}, sym_obs::VariableModel)
+ @assert sym_obs in m.g "$(sym_obs) is not observed."
+ @assert length(timeline) == length(observations)
+ nbr_observations = length(observations)
+
+ # Locations
+ locations = [:l0, :l1, :l2]
+
+ ## Invariant predicates
+ true_inv_predicate(x::Vector{Int}) = true
+ Λ_F = Dict(:l0 => true_inv_predicate, :l1 => true_inv_predicate,
+ :l2 => true_inv_predicate)
+
+ ## Init and final loc
+ locations_init = [:l0]
+ locations_final = [:l2]
+
+ map_var_automaton_idx = Dict{VariableAutomaton,Int}(:t => 1, :n => 2,
+ :d => 3, :idx => 4)
+ nbr_first_vars = 4
+ for i = 1:nbr_observations
+ map_var_automaton_idx[Symbol("tml_$(convert(Float64, i))")] = nbr_first_vars + i
+ end
+ for i = 1:nbr_observations
+ map_var_automaton_idx[Symbol("y_$(convert(Float64, i))")] = nbr_first_vars + nbr_observations + i
+ end
+
+ vector_flow = zeros(nbr_first_vars + 2*nbr_observations)
+ vector_flow[1] = 1.0
+ flow = Dict{Location,Vector{Float64}}(:l0 => vector_flow,
+ :l1 => vector_flow,
+ :l2 => vector_flow)
+
+ ## 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]
+ nbr_rand = rand(1:1000)
+ basename_func = "$(replace(m.name, ' '=>'_'))_$(nbr_rand)"
+ basename_func = replace(basename_func, '-'=>'_')
+
+ # l0 loc
+ # l0 => l1
+ sym_func_us_l0l1_1 = Symbol("us_eucl_dist_$(basename_func)_l0l1_1!")
+ str_us_l0l1_1 = "
+ @everywhere $(sym_func_us_l0l1_1)(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) = \n
+ (S.loc = :l1; \n
+ S[:n] = x[$(idx_obs_var)];\n
+ S[:d] = 0.0;\n
+ S[:idx] = 1.0)"
+ eval(Meta.parse(str_us_l0l1_1))
+ edge1 = Edge([nothing], getfield(Main, :cc_eucl_dist_aut_l0l1_1), getfield(Main, sym_func_us_l0l1_1))
+ map_edges[:l0][:l1] = [edge1]
+
+ # l1 loc
+ # l1 => l1
+ sym_func_us_l1l2_2 = Symbol("us_eucl_dist_$(basename_func)_l1l2_2!")
+ str_us_l1l2_2 = "
+ @everywhere $(sym_func_us_l1l2_2)(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =\n
+ (S[:n] = x[$(idx_obs_var)])"
+ eval(Meta.parse(str_us_l1l2_2))
+ edge1 = Edge([nothing], getfield(Main, :cc_eucl_dist_aut_l1l1_1), getfield(Main, :us_eucl_dist_aut_l1l1_1!))
+ edge2 = Edge([:ALL], getfield(Main, :cc_eucl_dist_aut_l1l1_2), getfield(Main, sym_func_us_l1l2_2))
+ map_edges[:l1][:l1] = [edge1, edge2]
+
+ # l1 => l2
+ edge1 = Edge([nothing], getfield(Main, :cc_eucl_dist_aut_l1l2_1), getfield(Main, :us_eucl_dist_aut_l1l2_1!))
+ map_edges[:l1][:l2] = [edge1]
+
+ ## Constants
+ constants = Dict{Symbol,Float64}(:nbr_obs => nbr_observations)
+ for i = 1:nbr_observations
+ constants[Symbol("tml_$(convert(Float64, i))")] = timeline[i]
+ constants[Symbol("y_$(convert(Float64, i))")] = observations[i]
+ end
+
+ 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
+
diff --git a/automata/period_automaton.jl b/automata/period_automaton.jl
new file mode 100644
index 0000000000000000000000000000000000000000..fe6d7f430e3d1089adb3cd2297f11a58f7b11c49
--- /dev/null
+++ b/automata/period_automaton.jl
@@ -0,0 +1,258 @@
+
+## Edge functions
+
+# l0 loc
+# * l0 => l0
+@everywhere cc_aut_per_l0l0_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+true
+@everywhere us_aut_per_l0l0_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(nothing)
+
+# * l0 => l0prime
+@everywhere cc_aut_per_l0l0prime_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:t] >= constants[:initT]
+@everywhere us_aut_per_l0l0prime_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S.loc = :l0prime)
+
+# * l0 => low
+@everywhere cc_aut_per_l0low_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:t] >= constants[:initT]
+@everywhere us_aut_per_l0low_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S.loc = :low;
+ S[:t] = 0.0;
+ S[:top] = 0.0;
+ S[:n] = -1)
+
+# l0prime
+# * l0prime => l0prime
+@everywhere cc_aut_per_l0primel0prime_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+true
+@everywhere us_aut_per_l0primel0prime_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(nothing)
+
+# * l0prime => low
+@everywhere cc_aut_per_l0primelow_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+true
+@everywhere us_aut_per_l0primelow_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S.loc = :low;
+ S[:t] = 0.0;
+ S[:top] = 0.0;
+ S[:n] = -1)
+
+# low
+# * low => low
+@everywhere cc_aut_per_lowlow_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:n] < constants[:N]
+@everywhere us_aut_per_lowlow_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(nothing)
+
+# * low => mid
+@everywhere cc_aut_per_lowmid_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:n] < constants[:N]
+@everywhere us_aut_per_lowmid_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S.loc = :mid)
+
+# * low => final
+@everywhere cc_aut_per_lowfinal_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:n] == constants[:N]
+@everywhere us_aut_per_lowfinal_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S.loc = :final)
+
+# mid
+# * mid => mid
+@everywhere cc_aut_per_midmid_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:n] < constants[:N]
+@everywhere us_aut_per_midmid_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(nothing)
+
+# * mid => low
+@everywhere cc_aut_per_midlow_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:n] < constants[:N] &&
+S[:top] == 0.0
+@everywhere us_aut_per_midlow_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S.loc = :low)
+
+@everywhere cc_aut_per_midlow_2(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:n] == -1.0 &&
+S[:top] == 1.0
+@everywhere us_aut_per_midlow_2!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S.loc = :low;
+ S[:n] += 1;
+ S[:t] = 0.0;
+ S[:top] = 0.0)
+
+@everywhere cc_aut_per_midlow_3(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(0 <= S[:n] <= 1) &&
+S[:top] == 1.0
+@everywhere us_aut_per_midlow_3!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S.loc = :low;
+ S[:n] += 1;
+ S[:top] = 0.0;
+ S[:mean_tp] = (S[:mean_tp] * (S[:n]-1) + S[:tp]) / S[:n];
+ S[:tp] = 0.0)
+
+@everywhere cc_aut_per_midlow_4(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(2 <= S[:n] < constants[:N]) &&
+S[:top] == 1.0
+@everywhere us_aut_per_midlow_4!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S.loc = :low;
+ S[:n] += 1;
+ S[:top] = 0.0;
+ S[:mean_tp] = (S[:mean_tp] * (S[:n]-1) + S[:tp]) / S[:n];
+ S[:var_tp] = (S[:var_tp] * (S[:n]-1) + (S[:mean_tp]-S[:tp])^2) / S[:n];
+ S[:tp] = 0.0)
+
+# * mid => high
+@everywhere cc_aut_per_midhigh_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:n] < constants[:N]
+@everywhere us_aut_per_midhigh_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S.loc = :high;
+ S[:top] = 1.0)
+
+# * mid => final
+@everywhere cc_aut_per_midfinal_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:n] == constants[:N]
+@everywhere us_aut_per_midfinal_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S.loc = :final)
+
+# high
+# * high => high
+@everywhere cc_aut_per_highhigh_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:n] < constants[:N]
+@everywhere us_aut_per_highhigh_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(nothing)
+
+# * high => mid
+@everywhere cc_aut_per_highmid_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:n] < constants[:N]
+@everywhere us_aut_per_highmid_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S.loc = :mid)
+
+# * high => final
+@everywhere cc_aut_per_highfinal_1(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+S[:n] == constants[:N]
+@everywhere us_aut_per_highfinal_1!(S::StateLHA, constants::Dict{Symbol,Float64}, x::Vector{Int}, p::Vector{Float64}) =
+(S.loc = :final)
+
+function create_period_automaton(m::ContinuousTimeModel, L::Float64, H::Float64, N::Int, sym_obs::VariableModel;
+ initT::Float64 = 0.0)
+ @assert sym_obs in m.g "$(sym_obs) is not observed."
+ N = convert(Float64, N)
+ nbr_rand = rand(1:1000)
+ basename_func = "$(replace(m.name, ' '=>'_'))_$(nbr_rand)"
+ basename_func = replace(basename_func, '-'=>'_')
+
+ ## Locations
+ locations = [:l0, :l0prime, :low, :mid, :high, :final]
+
+ ## Invariant predicates
+ idx_sym_obs = getfield(m, :map_var_idx)[sym_obs]
+ sym_name_L = Symbol("val_L_aut_per_$(basename_func)")
+ sym_name_H = Symbol("val_H_aut_per_$(basename_func)")
+
+ @everywhere true_predicate(x::Vector{Int}) = true
+ @everywhere low_predicate(x::Vector{Int}) = x[$(Meta.quot(idx_sym_obs))] <= $L
+ @everywhere not_low_predicate(x::Vector{Int}) = !low_predicate(x)
+ @everywhere mid_predicate(x::Vector{Int}) = $L < x[$(Meta.quot(idx_sym_obs))] < $H
+ @everywhere high_predicate(x::Vector{Int}) = x[$(Meta.quot(idx_sym_obs))] >= $H
+ #=
+ eval(Meta.parse("@everywhere true_predicate(x::Vector{Int}) = true"))
+ eval(Meta.parse("@everywhere low_predicate(x::Vector{Int}) = x[$(Meta.quot(idx_sym_obs))] <= $L"))
+ eval(Meta.parse("@everywhere not_low_predicate(x::Vector{Int}) = !low_predicate(x)"))
+ eval(Meta.parse("@everywhere mid_predicate(x::Vector{Int}) = L < x[$(Meta.quot(idx_sym_obs))] < $H"))
+ eval(Meta.parse("@everywhere high_predicate(x::Vector{Int}) = x[$(Meta.quot(idx_sym_obs))] >= $H"))
+ =#
+
+ Λ_F = Dict(:l0 => getfield(Main, :true_predicate), :l0prime => getfield(Main, :not_low_predicate),
+ :low => getfield(Main, :low_predicate), :mid => getfield(Main, :mid_predicate),
+ :high => getfield(Main, :high_predicate), :final => getfield(Main, :true_predicate))
+
+ ## Init and final loc
+ locations_init = [:l0]
+ locations_final = [:final]
+
+ ##Â Map of automaton variables
+ map_var_automaton_idx = Dict{VariableAutomaton,Int}(:t => 1, :n => 2, :top => 3, :tp => 4,
+ :mean_tp => 5, :var_tp => 6)
+
+ flow = Dict{Location,Vector{Float64}}(:l0 => [1.0,0.0,0.0,0.0,0.0,0.0],
+ :l0prime => [1.0,0.0,0.0,0.0,0.0,0.0],
+ :low => [1.0,0.0,0.0,1.0,0.0,0.0],
+ :mid => [1.0,0.0,0.0,1.0,0.0,0.0],
+ :high => [1.0,0.0,0.0,1.0,0.0,0.0],
+ :final => [1.0,0.0,0.0,0.0,0.0,0.0])
+
+ ## Edges
+ map_edges = Dict{Location, Dict{Location, Vector{Edge}}}()
+ for loc in locations
+ map_edges[loc] = Dict{Location, Vector{Edge}}()
+ end
+
+ # l0 loc
+ # * l0 => l0
+ edge_1 = Edge([:ALL], getfield(Main, :cc_aut_per_l0l0_1), getfield(Main, :us_aut_per_l0l0_1!))
+ map_edges[:l0][:l0] = [edge_1]
+ # * l0 => l0prime
+ edge_1 = Edge([nothing], getfield(Main, :cc_aut_per_l0l0prime_1), getfield(Main, :us_aut_per_l0l0prime_1!))
+ map_edges[:l0][:l0prime] = [edge_1]
+ # * l0 => low
+ edge_1 = Edge([nothing], getfield(Main, :cc_aut_per_l0low_1), getfield(Main, :us_aut_per_l0low_1!))
+ map_edges[:l0][:low] = [edge_1]
+
+ # l0prime
+ # * l0prime => l0prime
+ edge_1 = Edge([:ALL], getfield(Main, :cc_aut_per_l0primel0prime_1), getfield(Main, :us_aut_per_l0primel0prime_1!))
+ map_edges[:l0prime][:l0prime] = [edge_1]
+ # * l0prime => low
+ edge_1 = Edge([nothing], getfield(Main, :cc_aut_per_l0primelow_1), getfield(Main, :us_aut_per_l0primelow_1!))
+ map_edges[:l0prime][:low] = [edge_1]
+
+ # low
+ # * low => low
+ edge_1 = Edge([:ALL], getfield(Main, :cc_aut_per_lowlow_1), getfield(Main, :us_aut_per_lowlow_1!))
+ map_edges[:low][:low] = [edge_1]
+ # * low => mid
+ edge_1 = Edge([:ALL], getfield(Main, :cc_aut_per_lowmid_1), getfield(Main, :us_aut_per_lowmid_1!))
+ map_edges[:low][:mid] = [edge_1]
+ # * low => final
+ edge_1 = Edge([nothing], getfield(Main, :cc_aut_per_lowfinal_1), getfield(Main, :us_aut_per_lowfinal_1!))
+ map_edges[:low][:final] = [edge_1]
+
+ # mid
+ # * mid => mid
+ edge_1 = Edge([:ALL], getfield(Main, :cc_aut_per_midmid_1), getfield(Main, :us_aut_per_midmid_1!))
+ map_edges[:mid][:mid] = [edge_1]
+ # * mid => low
+ edge_1 = Edge([:ALL], getfield(Main, :cc_aut_per_midlow_1), getfield(Main, :us_aut_per_midlow_1!))
+ edge_2 = Edge([:ALL], getfield(Main, :cc_aut_per_midlow_2), getfield(Main, :us_aut_per_midlow_2!))
+ edge_3 = Edge([:ALL], getfield(Main, :cc_aut_per_midlow_3), getfield(Main, :us_aut_per_midlow_3!))
+ edge_4 = Edge([:ALL], getfield(Main, :cc_aut_per_midlow_4), getfield(Main, :us_aut_per_midlow_4!))
+ map_edges[:mid][:low] = [edge_1, edge_2, edge_3, edge_4]
+ # * mid => high
+ edge_1 = Edge([:ALL], getfield(Main, :cc_aut_per_midhigh_1), getfield(Main, :us_aut_per_midhigh_1!))
+ map_edges[:mid][:high] = [edge_1]
+ # * mid => final
+ edge_1 = Edge([nothing], getfield(Main, :cc_aut_per_midfinal_1), getfield(Main, :us_aut_per_midfinal_1!))
+ map_edges[:mid][:final] = [edge_1]
+
+ # high
+ # * high => high
+ edge_1 = Edge([:ALL], getfield(Main, :cc_aut_per_highhigh_1), getfield(Main, :us_aut_per_highhigh_1!))
+ map_edges[:high][:high] = [edge_1]
+ # * high => mid
+ edge_1 = Edge([:ALL], getfield(Main, :cc_aut_per_highmid_1), getfield(Main, :us_aut_per_highmid_1!))
+ map_edges[:high][:mid] = [edge_1]
+ # * high => final
+ edge_1 = Edge([nothing], getfield(Main, :cc_aut_per_highfinal_1), getfield(Main, :us_aut_per_highfinal_1!))
+ map_edges[:high][:final] = [edge_1]
+
+ ## Constants
+ constants = Dict{Symbol,Float64}(:N => N, :L => L, :H => H, :initT => initT)
+
+ A = LHA("Period", 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_period_automaton
+
diff --git a/examples/scripts/doping_3way_oscillator.jl b/examples/scripts/doping_3way_oscillator.jl
new file mode 100644
index 0000000000000000000000000000000000000000..cbfba14af7fad35a5e19e77e5207969620e3c996
--- /dev/null
+++ b/examples/scripts/doping_3way_oscillator.jl
@@ -0,0 +1,16 @@
+
+using MarkovProcesses
+load_plots()
+
+load_model("doping_3way_oscillator")
+load_automaton("period_automaton")
+A_per = create_period_automaton(doping_3way_oscillator, 300.0, 360.0, 5, :A)
+
+sync_doping = doping_3way_oscillator * A_per
+set_time_bound!(sync_doping, 0.1)
+set_x0!(doping_3way_oscillator, [:A, :B, :C, :DA, :DB, :DC], [333, 333, 333, 10, 10, 10])
+σ = simulate(sync_doping)
+plot(σ; A = A_per, filename = "traj_full.png")
+plot(σ, :A; A = A_per, filename = "traj_A.png")
+plot_periodic_trajectory(A_per, σ, :A, filename = "traj_automaton.png")
+