ABMEv_plot.jl 7.74 KB
Newer Older
1
using RecipesBase
2
3
using Colors
import KernelDensity:kde,pdf
Victor's avatar
Victor committed
4
5
6
7
8
"""
    function plot(world::Array{U},p;what=["x","H"],trait = 1,tplot = false) where U <: Union{Missing,Agent}

# ARGS
- `what = ["x","H"]`: the plots you want to obtain
9
- `trait = 1`: the trait that will plotted regarding what you asked. `trait = 0` will plot the geotrait
10
- `tplot = 0` used when calling xs, as it plots a snapshot of the world at a particular time
Victor's avatar
Victor committed
11
It should correspond to an integer, as it indexes the column to plot
12
13
14
15

# Options available
- `"x"`
- `"xs"`
Victor's avatar
Victor committed
16
17
18
"""

@recipe function plot(world::Array{U},p;what=["x","H"],trait = 1,tplot = 0) where U <: Union{Missing,Agent}
19
    tot_dim = size(world,2)*size(world,1)
20
    # We reduce time interval if it is too big
Victor's avatar
Victor committed
21
22
23
24
25
26
    # if tot_dim > 1e6 && size(world,2) >= 200
    #     p = copy(p)
    #     idx_reduced = floor.(Int,range(1,size(world,2),length = 200))
    #     p["tspan" ] = p["tspan"][idx_reduced]
    #     world = world[:,idx_reduced]
    # end
27
28
29
    # second condition is to make sure that the world corresponds to all the time steps
    # If not, then we can not get "x"
    if count(ismissing,world) > 0 && length(p["tspan"]) == size(world,2)
30
31
32
33
34
        tspan_ar = vcat([p["tspan"][i]*ones(Int(p["NMax"] - count(ismissing,world[:,i]))) for i in 1:length(p["tspan"]) ]...);
    else
        tspan_ar = repeat(p["tspan"],inner = size(world,1))
    end
    # tspan = Float64.(tspan)
35
    tend = p["tspan"][tplot > 0 ? tplot : length(p["tspan"])]
36
    world_sm = clean_world(world)
37
    if "x" in what
38
39
        d_i = []
        for i in 1:size(world,2)
40
            x = get_x(clean_world(world[:,i]),p["tspan"][i],trait)[:]
41
42
            append!(d_i,pdf(kde(x),x))
        end
43
        @series begin
44
45
46
            if length(world_sm) !== length(tspan_ar[:])
                throw(DimensionMismatch("You want to plot a world with missing data"))
            end
47
        xarray = get_x.(world_sm,tspan_ar[:],trait)
48
            seriestype := :scatter
49
            markercolor := eth_grad_small[d_i ./ maximum(d_i)]
50
            markerstrokewidth := 0
Victor's avatar
Victor committed
51
            seriesalpha :=1.
Victor's avatar
Victor committed
52
53
            # xlabel := "time"
            # ylabel := "trait value"
54
            label := ""
55
56
            grid := false
            # markersize := 2.3/1000*size(world_sm,1)
57
            tspan_ar[:],xarray[:]
58
59
        end
    end
60
61
62
    # we use this for discrete agents
    # world should be a one dimensional vector, corresponding to one time step only
    if "xs" in what
Victor's avatar
Victor committed
63
        d_i = []; xt_array = []; x1_array = []
64
65
        world_df_all = world2df(clean_world(world[:, tplot > 0 ? tplot : size(world,2) ]),tend,true)
        world_df_g = groupby(world_df_all,:x1)
66
        for world_df in world_df_g
67
            if trait == 0
68
                x = Float64.(world_df.g)
69
70
71
72
73
            else
                # fitness occupies first spot
                x = world_df[:,trait+1] ;
            end
            x1 =  world_df.x1;
74
            append!(d_i,pdf(kde(x),x))
Victor's avatar
Victor committed
75
76
            append!(xt_array,x)
            append!(x1_array,x1)
77
78
79
80
81
82
83
        end
        # TODO: we stopped here
        @series begin
            seriestype := :scatter
            markercolor := eth_grad_small[d_i ./ maximum(d_i)]
            # markercolor := :blue
            markerstrokewidth := 0
84
            # seriesalpha := 1.
Victor's avatar
Victor committed
85
            xaxis := "geographical position"
86
            xticks :=  sort!(unique(world_df_all.x1))
Victor's avatar
Victor committed
87
            yaxis := "trait value"
88
89
            label := ""
            grid := false
Victor's avatar
Victor committed
90
            # marker := (:rect,20,1.)
Victor's avatar
Victor committed
91
            x1_array[:],xt_array[:]
92
93
        end
    end
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
    if "gs" in what
        _world = clean_world(world[:, tplot > 0 ? tplot : length(p["tspan"]) ])
        y = get_x(_world,tend,2)[:]
        x = get_x(_world,tend,0)[:]
        X = hcat(x,y)
        d = kde(X)
        # by density
        d_i = diag(pdf(d,X[:,1],X[:,2]))
        # by value
        # d_i = y
        d_i = (d_i .- minimum(d_i)) ./ (maximum(d_i) .- minimum(d_i))
        # TODO: we stopped here
        @series begin
            seriestype := :scatter
            markercolor := eth_grad_small[d_i]
            # markercolor := :blue
            markerstrokewidth := 0
            # seriesalpha := 1.
            xaxis := "geotrait"
            yaxis := "trait value"
            label := ""
            grid := false
            # marker := (:rect,20,1.)
            x,y
        end
    end
120
    if "3dgeo" in what
121
122
123
124
125
126
127
128
129
130
131
132
133
        d_i = []
        for i in 1:size(world,2)
            _world = clean_world(world[:,i])
            x = get_x(_world,p["tspan"][i],2)[:]
            y = get_x(_world,p["tspan"][i],0)[:]
            X = hcat(x,y)
            # d = kde(X)
            # di_temp = diag(pdf(d,X[:,1],X[:,2]))
            di_temp = y
            di_temp = (di_temp .- minimum(di_temp)) ./ (maximum(di_temp) .- minimum(di_temp))
            # here we normalise with respect to maximum value at each time step
            append!(d_i,di_temp)
        end
134
        @series begin
135
136
        xarray = get_geo.(world_sm,tspan_ar)
        yarray = get_x(world_sm,2)
137
            seriestype := :scatter3d
138
            markercolor := eth_grad_std[d_i ./ 1.]
139
            markerstrokewidth := 0
Victor's avatar
Victor committed
140
            seriesalpha :=.1
141
142
143
144
            xlabel := "time"
            ylabel := "geotrait"
            zlabel := "trait value"
            label := ""
145
            # markersize := 2.3/1000*size(world_sm,1)
146
            tspan_ar,xarray[:],yarray[:]
147
148
149
        end
    end
    if "3d" in what
150
151
152
153
154
155
156
157
158
159
        d_i = []
        for i in 1:size(world,2)
            x = get_x(clean_world(world[:,i]),p["tspan"][i],1)[:]
            y = get_x(clean_world(world[:,i]),p["tspan"][i],2)[:]
            X = hcat(x,y)
            d = kde(X)
            di_temp = diag(pdf(d,X[:,1],X[:,2]))
            di_temp = (di_temp .- minimum(di_temp)) ./ (maximum(di_temp) .- minimum(di_temp))
            append!(d_i,di_temp)
        end
160
        @series begin
161
        xarray = get_x(world_sm,1)
162
        yarray = get_x(world_sm,2)
163
            seriestype := :scatter3d
164
            markercolor := eth_grad_small[d_i ./ maximum(d_i)]
165
            markerstrokewidth := 0
Victor's avatar
Victor committed
166
            seriesalpha :=.1
167
168
169
170
            xlabel := "time"
            ylabel := "position"
            zlabel := "trait value"
            label := ""
171
            # markersize := 2.3/1000*size(world_sm,1)
172
            tspan_ar,xarray[:],yarray[:]
173
174
175
176
        end
    end
    # if "H" in what
    #     @series begin
177
    #         x = get_x.(world_sm,trait)
178
179
180
181
182
183
184
185
186
187
188
189
190
    #         linewidth := 2
    #         seriestype := :line
    #         label := "Interconnectedness"
    #         tspan,N^2 / 2 .* [H_discrete(x[:,i]) for i in tspan]
    #     end
    # end
    if "var" in what
        @series begin
            linewidth := 2
            seriestype := :line
            label := "Variance"
            xlabel := "Time"
            ylabel := "Variance"
191
            p["tspan"],var(world_sm,trait=trait)[:]
192
193
194
195
196
197
198
199
200
        end
    end
    if "vargeo" in what
        @series begin
            linewidth := 2
            seriestype := :line
            label := "Variance of geotrait"
            xlabel := "Time"
            ylabel := "Variance"
201
            p["tspan"],i->first(covgeo(world_sm[:,Int(i)]))
202
203
        end
    end
204
205
206
207
208
209
210
    # if "density_t" in what
    #     @series begin
    #         linewidth := 2
    #         seriestype := :plot3d
    #         label := "Variance of geotrait"
    #         xlabel := "Time"
    #         ylabel := "Variance"
211
    #         p["tspan"],i->first(covgeo(world_sm[:,Int(i)]))
212
213
    #     end
    # end
214
end
Victor's avatar
Victor committed
215
216
217
218
219
220
221


import Plots:cgrad
# asymmetry towards red, blue is only a fifth of the color range
const eth_grad_small = cgrad([colorant"#1F407A", RGB(0.671,0.851,0.914),RGB(1.0,1.0,0.749), RGB(0.992,0.682,0.38),RGB(0.647,0.0,0.149),],[.0,.1])
# symmetry between red and blue
const eth_grad_std = cgrad([colorant"#1F407A", RGB(0.671,0.851,0.914),RGB(1.0,1.0,0.749), RGB(0.992,0.682,0.38),RGB(0.647,0.0,0.149),],[.0,1.])