plots.jl 5.76 KB
Newer Older
1

2
import Plots: plot, plot!, scatter!, hline!, Shape, text
3
import Plots: current, palette, display, png, close, savefig
4

5
6
7
8
9
"""
    `plot(σ, var...; plot_transitions=false)`

Plot a simulated trajectory σ. var... is a tuple of stirng variables.
`plot(σ)` will plot all the variables simulated in σ 
10
whereas `plot(σ, :I, :R)` only plots the variables I and R of the trajectory (if it exists).
11
12
13
If `plot_transitions=true`, a marker that corresponds to a transition of the model will be plotted
at each break of the trajectory.
"""
14
function plot(σ::AbstractTrajectory, vars::VariableModel...; plot_transitions::Bool = false, A::Union{Nothing,LHA} = nothing, filename::String = "")
15
16
17
18
19
20
21
22
23
24
    # Setup 
    palette_tr = palette(:default)
    l_tr = unique(transitions(σ))
    map_tr_color(tr) = palette_tr[findfirst(x->x==tr, l_tr)]
    to_plot = vars
    if length(vars) ==  0
        to_plot = get_obs_var(σ)
    end
    
    # Plots
25
    p = plot(title = "Trajectory of $(σ.m.name) model", palette = :lightrainbow, legend = :outertopright, background_color_legend=:transparent, dpi = 480)
26
    for var in to_plot
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
27
        @assert var in get_obs_var(σ) "Variable $var is not observed." 
28
29
        plot!(p, times(σ), σ[var], 
              xlabel = "Time", ylabel = "Number of species",
30
              label = "$var",
31
32
33
34
35
36
              linetype=:steppost)
    end
    if plot_transitions
        for (i, var) in enumerate(to_plot)
            for tr in l_tr
                idx_tr = findall(x->x==tr, transitions(σ))
37
                label = (tr == nothing || i > 1) ? "" : "$tr"
38
39
40
41
42
43
44
45
                alpha = (tr == nothing) ? 0.0 : 0.5
                scatter!(p, times(σ)[idx_tr], σ[var][idx_tr], label=label, 
                         markershape=:cross, markeralpha=alpha, 
                         markersize = 2,
                         markercolor=palette_tr[findfirst(x->x==tr, l_tr)])
            end
        end
    end
46
47
48
    if A != nothing
        plot!(A)
    end
49
50
51
    if filename == ""
        display(p)
    else
52
        savefig(p, filename)
53
54
55
    end
end

Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
56
57
function plot!(A::LHA; label::String = "")
    label_region = (label == "") ? "Region LHA" : label
58
    if A.name in ["F property", "G property"]
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
59
60
        if haskey(A.constants, :x1) && haskey(A.constants, :x2) && haskey(A.constants, :t1) && haskey(A.constants, :t2) 
            x1, x2, t1, t2 = A.constants[:x1], A.constants[:x2], A.constants[:t1], A.constants[:t2] 
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
61
            plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5, label = label_region)
62
63
64
        end
    end
    if A.name in ["G and F property"]    
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
65
66
        if haskey(A.constants, :x1) && haskey(A.constants, :x2) && haskey(A.constants, :t1) && haskey(A.constants, :t2) 
            x1, x2, t1, t2 = A.constants[:x1], A.constants[:x2], A.constants[:t1], A.constants[:t2] 
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
67
            plot!(Shape([(t1,x1), (t1,x2), (t2,x2), (t2,x1), (t1,x1)]), opacity = 0.5, label = label_region)
68
        end
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
69
70
        if haskey(A.constants, :x3) && haskey(A.constants, :x4) && haskey(A.constants, :t3) && haskey(A.constants, :t4) 
            x3, x4, t3, t4 = A.constants[:x3], A.constants[:x4], A.constants[:t3], A.constants[:t4] 
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
71
            plot!(Shape([(t3,x3), (t3,x4), (t4,x4), (t4,x3), (t3,x3)]), opacity = 0.5, label = label_region)
72
73
        end
    end
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
74
    label_region = (label == "") ? "L, H" : label
75
    if A.name == "Period"
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
76
        hline!([A.constants[:L], A.constants[:H]], label = label_region, linestyle = :dot)
77
78
79
80
    end
end

# For tests purposes
81
82
function plot_periodic_trajectory(A::LHA, σ::SynchronizedTrajectory, sym_obs::Symbol; verbose = false, filename::String = "")
    @assert sym_obs in get_obs_var(σ) "Variable is not observed in the model"
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
    @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)
    nbr_states = length_states(σ)
    locations_trajectory = Vector{Location}(undef, nbr_states)
    locations_trajectory[1] = Sn.loc
    idx_n = [1]
    values_n = [Sn[:n]]
    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]
            push!(idx_n, n)
            push!(values_n, Sn[:n])
        end
        if Snplus1.loc in A_new.locations_final 
            break 
        end
108
    end
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
109
    p = plot(title = "Oscillatory trajectory of $(σ.m.name) model", palette = :lightrainbow,
110
111
112
113
114
115
             background_color_legend=:transparent, dpi = 480, legend = :outertopright) #legendfontsize, legend
    colors_loc = Dict(:l0 => :silver, :l0prime => :silver, :final => :black,
                      :low => :skyblue2, :mid => :orange, :high => :red)
    loc_to_color(loc::Symbol) = colors_loc[loc]
    for loc in A.locations
        idx_locations = findall(x->x==loc, locations_trajectory)
116
        label_state = (loc in [:low, :mid, :high]) ? "$sym_obs ($loc loc)" : ""
117
118
        scatter!(p, times(σ)[idx_locations], σ[sym_obs][idx_locations], color = colors_loc[loc], 
                    markersize = 1.0, markershape = :cross, 
119
                    label = label_state, xlabel = "Time", ylabel = "Species $sym_obs")
120
    end
121
122
123
124
    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,
                             markershape = :utriangle, markersize = 3, label = "n")
    hline!(p, [A.constants[:L], A.constants[:H]], label = "L, H", color = :grey; linestyle = :dot)
125
126
127
128
    
    if filename == ""
        display(p)
    else
129
        savefig(p, filename)
130
    end
131
132
end

133
export plot, plot!, plot_periodic_trajectory
134