ABMEv_plot.jl 7.21 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
        xarray = get_x.(world_sm,tspan_ar[:],trait)
45
            seriestype := :scatter
46
47
            markercolor := eth_grad_small[d_i ./ maximum(d_i)]
            # markercolor := :blue
48
            markerstrokewidth := 0
Victor's avatar
Victor committed
49
            seriesalpha :=1.
50
51
52
            xlabel := "time"
            ylabel := "trait value"
            label := ""
53
54
            grid := false
            # markersize := 2.3/1000*size(world_sm,1)
55
            tspan_ar[:],xarray[:]
56
57
        end
    end
58
59
60
    # 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
61
        d_i = []; xt_array = []; x1_array = []
62
63
        world_df_all = world2df(clean_world(world[:, tplot > 0 ? tplot : size(world,2) ]),tend,true)
        world_df_g = groupby(world_df_all,:x1)
64
        for world_df in world_df_g
65
            if trait == 0
66
                x = Float64.(world_df.g)
67
68
69
70
71
            else
                # fitness occupies first spot
                x = world_df[:,trait+1] ;
            end
            x1 =  world_df.x1;
72
            append!(d_i,pdf(kde(x),x))
Victor's avatar
Victor committed
73
74
            append!(xt_array,x)
            append!(x1_array,x1)
75
76
77
78
79
80
81
        end
        # TODO: we stopped here
        @series begin
            seriestype := :scatter
            markercolor := eth_grad_small[d_i ./ maximum(d_i)]
            # markercolor := :blue
            markerstrokewidth := 0
82
            # seriesalpha := 1.
Victor's avatar
Victor committed
83
            xaxis := "geographical position"
84
            xticks :=  sort!(unique(world_df_all.x1))
Victor's avatar
Victor committed
85
            yaxis := "trait value"
86
87
            label := ""
            grid := false
88
            marker := (:rect,20,1.)
Victor's avatar
Victor committed
89
            x1_array[:],xt_array[:]
90
91
        end
    end
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
    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
118
    if "3dgeo" in what
119
120
121
122
123
124
125
126
127
128
129
130
131
        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
132
        @series begin
133
134
        xarray = get_geo.(world_sm,tspan_ar)
        yarray = get_x(world_sm,2)
135
            seriestype := :scatter3d
136
            markercolor := eth_grad_std[d_i ./ 1.]
137
            markerstrokewidth := 0
Victor's avatar
Victor committed
138
            seriesalpha :=.1
139
140
141
142
            xlabel := "time"
            ylabel := "geotrait"
            zlabel := "trait value"
            label := ""
143
            # markersize := 2.3/1000*size(world_sm,1)
144
            tspan_ar,xarray[:],yarray[:]
145
146
147
        end
    end
    if "3d" in what
148
149
150
151
152
153
154
155
156
157
        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
158
        @series begin
159
        xarray = get_x(world_sm,1)
160
        yarray = get_x(world_sm,2)
161
            seriestype := :scatter3d
162
            markercolor := eth_grad_small[d_i ./ maximum(d_i)]
163
            markerstrokewidth := 0
Victor's avatar
Victor committed
164
            seriesalpha :=.1
165
166
167
168
            xlabel := "time"
            ylabel := "position"
            zlabel := "trait value"
            label := ""
169
            # markersize := 2.3/1000*size(world_sm,1)
170
            tspan_ar,xarray[:],yarray[:]
171
172
173
174
        end
    end
    # if "H" in what
    #     @series begin
175
    #         x = get_x.(world_sm,trait)
176
177
178
179
180
181
182
183
184
185
186
187
188
    #         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"
189
            p["tspan"],var(world_sm,trait=trait)[:]
190
191
192
193
194
195
196
197
198
        end
    end
    if "vargeo" in what
        @series begin
            linewidth := 2
            seriestype := :line
            label := "Variance of geotrait"
            xlabel := "Time"
            ylabel := "Variance"
199
            p["tspan"],i->first(covgeo(world_sm[:,Int(i)]))
200
201
        end
    end
202
203
204
205
206
207
208
    # if "density_t" in what
    #     @series begin
    #         linewidth := 2
    #         seriestype := :plot3d
    #         label := "Variance of geotrait"
    #         xlabel := "Time"
    #         ylabel := "Variance"
209
    #         p["tspan"],i->first(covgeo(world_sm[:,Int(i)]))
210
211
    #     end
    # end
212
end