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