diff --git a/automata/period_automaton.jl b/automata/period_automaton.jl index cfd490cdbc103a3a33bc8a544089cf2ecc442ea2..89c6ee276d57312a9a8b33c19cfd8ae8811b1174 100644 --- a/automata/period_automaton.jl +++ b/automata/period_automaton.jl @@ -153,8 +153,8 @@ function create_period_automaton(m::ContinuousTimeModel, L::Float64, H::Float64, @everywhere $(func_name(:us, :mid, :low, 2))(S::StateLHA, x::Vector{Int}, p::Vector{Float64}) = (setfield!(S, :loc, Symbol("low")); setindex!(getfield(S, :values), get_value(S, $(get_idx_var(:n))) + 1, $(get_idx_var(:n))); - setindex!(getfield(S, :values), 0.0, $(get_idx_var(:t))); - setindex!(getfield(S, :values), 0.0, $(get_idx_var(:top)))) + setindex!(getfield(S, :values), 0.0, $(get_idx_var(:top))); + setindex!(getfield(S, :values), 0.0, $(get_idx_var(:tp)))) @everywhere $(func_name(:cc, :mid, :low, 3))(S::StateLHA, x::Vector{Int}, p::Vector{Float64}) = (0 <= get_value(S, $(get_idx_var(:n))) <= 1) && @@ -181,7 +181,8 @@ function create_period_automaton(m::ContinuousTimeModel, L::Float64, H::Float64, setindex!(getfield(S, :values), g_var_tp(get_value(S, $(get_idx_var(:var_tp))), get_value(S, $(get_idx_var(:mean_tp))), get_value(S, $(get_idx_var(:tp))), - get_value(S, $(get_idx_var(:n)))), $(get_idx_var(:var_tp)))) + get_value(S, $(get_idx_var(:n)))), $(get_idx_var(:var_tp))); + setindex!(getfield(S, :values), 0.0, $(get_idx_var(:tp)))) # * mid => high @everywhere $(func_name(:cc, :mid, :high, 1))(S::StateLHA, x::Vector{Int}, p::Vector{Float64}) = diff --git a/core/model.jl b/core/model.jl index 9a9922552953dde01ca427bb4af8dcc19576980f..4f33d5c275c0f43027b39b4d0fd96d1665f930b1 100644 --- a/core/model.jl +++ b/core/model.jl @@ -160,7 +160,7 @@ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Fl vec_x = zeros(Int, getfield(m, :dim_state)) l_t = Float64[0.0] l_tr = Transition[nothing] - Snplus1 = copy(Sn) + Snplus1 = deepcopy(Sn) # If x0 is absorbing if isabsorbing || isacceptedLHA _resize_trajectory!(full_values, times, transitions, 1) @@ -288,7 +288,7 @@ function volatile_simulate(product::SynchronizedModel; vec_x = zeros(Int, getfield(m, :dim_state)) l_t = Float64[0.0] l_tr = Transition[nothing] - Snplus1 = copy(Sn) + Snplus1 = deepcopy(Sn) # If x0 is absorbing if isabsorbing || isacceptedLHA if !isacceptedLHA diff --git a/core/plots.jl b/core/plots.jl index bfdf9fe8e1779400b344ba12517de2c1ea53333c..e23ba88ce1e7639138856786c57059e22c023f75 100644 --- a/core/plots.jl +++ b/core/plots.jl @@ -78,31 +78,32 @@ function plot!(A::LHA; label::String = "") end # For tests purposes -function plot_periodic_trajectory(A::LHA, σ::SynchronizedTrajectory, sym_obs::Symbol; verbose = false, filename::String = "") +function plot_periodic_trajectory(A::LHA, σ::SynchronizedTrajectory, sym_obs::Symbol; verbose = false, annot_size::Float64 = 6.0, filename::String = "") @assert sym_obs in get_obs_var(σ) "Variable is not observed in the model" @assert A.name in ["Period"] - A_new = LHA(A, (σ.m)._map_obs_var_idx) p_sim = (σ.m).p l_t = times(σ) l_tr = transitions(σ) - Sn = init_state(A_new, σ[1], l_t[1]) - Snplus1 = copy(Sn) + Sn = init_state(A, σ[1], l_t[1]) + Snplus1 = deepcopy(Sn) nbr_states = length_states(σ) locations_trajectory = Vector{Location}(undef, nbr_states) locations_trajectory[1] = Sn.loc idx_n = [1] values_n = [Sn[:n]] + values_tp = [Sn[:tp]] if verbose println("Init: ") end if verbose @show Sn end for n in 2:nbr_states - next_state!(Snplus1, A_new, σ[n], l_t[n], l_tr[n], Sn, σ[n-1], p_sim; verbose = verbose) - copyto!(Sn, Snplus1) - locations_trajectory[n] = Sn.loc - if Sn[:n] != values_n[end] + next_state!(Snplus1, A, σ[n], l_t[n], l_tr[n], Sn, σ[n-1], p_sim; verbose = verbose) + if Snplus1[:n] != values_n[end] push!(idx_n, n) - push!(values_n, Sn[:n]) + push!(values_n, Snplus1[:n]) + push!(values_tp, Sn[:tp]) end - if Snplus1.loc in A_new.locations_final + copyto!(Sn, Snplus1) + locations_trajectory[n] = Sn.loc + if Snplus1.loc in A.locations_final break end end @@ -118,8 +119,10 @@ function plot_periodic_trajectory(A::LHA, σ::SynchronizedTrajectory, sym_obs::S markersize = 1.0, markershape = :cross, label = label_state, xlabel = "Time", ylabel = "Species $sym_obs") end - annot_n = [(times(σ)[idx_n[i]], σ[sym_obs][idx_n[i]] - 10, text("n = $(values_n[i])", 6, :bottom)) for i = eachindex(idx_n)] - scatter!(p, times(σ)[idx_n], σ[sym_obs][idx_n], annotations = annot_n, + annot_n = [(times(σ)[idx_n[i]], σ[sym_obs][idx_n[i]] - 10, text("n = $(values_n[i])", annot_size, :top)) for i = eachindex(idx_n)] + annot_tp = [(times(σ)[idx_n[i]], σ[sym_obs][idx_n[i]] - 10, text("tp = $(round(values_tp[i], digits = 2))", annot_size, :bottom)) for i = eachindex(idx_n)] + annots = vcat(annot_n, annot_tp) + scatter!(p, times(σ)[idx_n], σ[sym_obs][idx_n], annotations = annots, markershape = :utriangle, markersize = 3, label = "n") hline!(p, [A.constants[:L], A.constants[:H]], label = "L, H", color = :grey; linestyle = :dot)