ABMEv_metrics.jl 7.01 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
import StatsBase:Histogram,fit,step,normalize

d1(x,y) = norm(x .- y,1)
d2(x,y) = norm(x .- y,2)
hamming(k::Array, l::Array) = sum(k .!= l)

"""
    H_discrete(s)
Interconnectedness measure as in Nordbotten 2018 for discrete setup
"""
function H_discrete(s)
    N = length(s)
    H = 0
    for x in s
        for y in s
            H += d2(x,y)
        end
    end
    return H / N^2
end


"""
    findclusters(v::Vector,allextrema =true)
Returns a tuple with the cluster mean and its associated weight
# Arguments

"""
function findclusters(v::Vector,allextrema =true)
    # this function fits a histogram to some distribution and extracts its local maxima
    h = fit(Histogram,v)
    s = step(h.edges...)
    x = normalize(h.weights,1)
    dx = x[2:end] .- x[1:end-1]
    if allextrema
        sdx = dx[2:end] .* dx[1:end-1]
        idx = findall(i -> i < 0, sdx)
    else
        idx = findall(i -> dx[i] > 0 && dx[i+1] < 0, 1:length(dx))
    end
    return collect(h.edges...)[idx .+ 1] .+ s/2, x[idx .+ 1]
end

import Statistics.var
45
# TODO: rename this to gamma diversity
46
47
48
49
50
51
52
53
"""
    var(world::Array{Agent};trait=:)
If trait = 0, returns the variance of the geotrait,
knowing that by default it is associated with position trait 1.
If trait > 0, returns the covariance matrix, with first row geotrait and second row trait.
# Arguments

"""
54
function var(world::World;trait=1)
55
    xarray = get_x(world,trait)
56
    return var(xarray,dims=1,corrected=false)
57
58
59
60
61
62
63
64
65
end
"""
    covgeo(world::Array{Agent,1},trait = 0)
If trait = 0, returns the variance of the geotrait,
knowing that by default it is associated with position trait 1.
If trait > 0, returns the covariance matrix, with first row geotrait and second row trait
# Arguments

"""
Victor's avatar
Victor committed
66
67
function covgeo(world::World,trait = 0)
    xarray = Float64.(get_geo(world))
68
    if trait > 0
Victor's avatar
Victor committed
69
        xstd = reshape(Float64.(get_x(world,trait)),size(world,1),size(world,2))
70
71
        xarray = hcat(xarray,xstd)
    end
72
    return cov(xarray,corrected=false)
73
74
75
76
77
78
79
80
end

"""
    function hamming(world::Array{Agent,1})
Returns a matrix H where H_ij = hamming(a_i,a_j).
The hamming distance is taken through the whole history
and functional space of the agents.
"""
Victor's avatar
Victor committed
81
82
83
84
function hamming(world::World) where T <: Int
    n = size(world)
    H = zeros(n,n)
    for (i,a) in enumerate(agents(world))
85
86
87
88
89
90
            for (j,b) in enumerate(world)
                    H[i,j] = hamming(get_xhist(a),get_xhist(b))
            end
    end
    return H
end
91
"""
92
    get_alpha_div(world::Array{U,1},t::Number,trait=1) where U <: Union{Missing,Agent}
Victor's avatar
Victor committed
93
Mean of the local variance of `trait` per patch. If trait=0, we get the geotrait
94
95
# Arguments
"""
Victor's avatar
Victor committed
96
97
98
function get_alpha_div(world::World,trait=1)
    xarray = get_xarray(world,true)
    g = groupby(x->x[1],collect(eachrow(xarray)))
Victor's avatar
Victor committed
99
    if trait == 0
Victor's avatar
Victor committed
100
        return mean([var(Float64.([x[end] for x in xp]),corrected=false) for xp in values(g)])
Victor's avatar
Victor committed
101
    else
Victor's avatar
Victor committed
102
        return mean([var(Float64.([x[trait] for x in xp]),corrected=false) for xp in values(g)])
Victor's avatar
Victor committed
103
    end
104
105
106
end

"""
107
    get_beta_div(world::Array{U,1},t::Number,trait=1) where U <: Union{Missing,Agent}
Victor's avatar
Victor committed
108
Variance of the mean of `trait` per patch
109
110
# Arguments
"""
Victor's avatar
Victor committed
111
112
113
function get_beta_div(world::World,trait=1)
    xarray = get_xarray(world,true)
    g = groupby(x->x[1],collect(eachrow(xarray)))
Victor's avatar
Victor committed
114
    if trait == 0
115
        # need to convert to Float64, otherwise infinite variance
Victor's avatar
Victor committed
116
        sbar_i = [mean(Float64.([x[end] for x in xp])) for xp in values(g)]
Victor's avatar
Victor committed
117
    else
Victor's avatar
Victor committed
118
        sbar_i = [mean(Float64.([x[trait] for x in xp])) for xp in values(g)]
Victor's avatar
Victor committed
119
    end
Victor's avatar
bug fix    
Victor committed
120
    return var(sbar_i,corrected=false)
121
end
122
123

"""
124
    function get_xhist_mat(world,trait=1,time = 0)
125
returns `xhist,ttot`, where `xhist` is a matrix with dimension `lenght(world)` x `length(thist)+1`,
126
127
128
which consists of geographical history of ancestors at every time step.
If `time` is specified and is higher that the highest time step in world,
then the last column of xhist corresponds to actual position of agents
129
"""
130
131
132
133
134
135

function get_xhist_mat(world,trait=1,time = 0)
        thist = vcat(get_thist.(world)...)
        ttot = sort!(unique(thist))
        xhist = zeros(length(world),length(ttot))
        for (i,a) in enumerate(world)
136
137
138
139
140
141
142
143
144
145
            _thist = get_thist(a)
            # Here we check that there is no redundancy of thist because of casting errors
            # In the future, we should remove this check, as the type of time has been set to Float64
            # Hence there should be no more problems of this type
            ttrue = _thist[2:end] .- _thist[1:end-1] .> 0
            if count(ttrue) < length(_thist) - 1
                _tt = vcat(ttrue,true)
                _thist = _thist[_tt] # we drop the position that was occupied for dt = .0
                _xhist = get_xhist(a,trait)[_tt]
            else
146
                _xhist = get_xhist(a,trait)
147
148
149
150
151
152
153
154
155
156
157
158
159
160
            end
            k = 0
            _l = length(_thist)
            for j in 1:(_l - 1)
                    xhist[i,j+k] = _xhist[j]
                            while _thist[j+1] > ttot[j+k+1]
                                    k+=1
                                    xhist[i,j+k] = _xhist[j]
                                    if j+k == length(ttot)
                                            break
                                    end
                            end
            end
            xhist[i,_l+k:end] .= _xhist[end]
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
        end
        tt = 0
        time > ttot[end] ? ttot = vcat(ttot,time) : tt = 1
        return xhist[:,1:end-tt],ttot
end

"""
    function get_dist_hist(a1,a2,dist,trait=1,time = 0)

Returns the integral of the distance `dist` through time of `trait` between `a1.x` and `a2.x`.

```math
\\int d(x_1,x_2,t)dt
```
"""
function get_dist_hist(a1,a2,dist,trait=1,time = 0)
177
178
        xhist,ttot = get_xhist_mat([a1,a2],trait,time)
        return sum(dist.(xhist[1,:],xhist[2,:]) .* (ttot[2:end] .- ttot[1:end-1]))
179
180
end

181
182
183
184
"""
    function truncvar(X::AbstractArray)
Returns the truncated variance of `X`.
"""
185

186
187
188
189
190
function truncvar(X::AbstractArray)
        N = length(X)
        c = [count(y->y!=x,X) for x  in unique(X)]
        return sum(c .* (N .-c)) /(2*N^2)
end
191

192
"""
193
194
195
    function get_pairwise_average_isolation(world;trait=1,trunc=false)
Returns the integrated pairwise squared distance between all agents of `world` wrt `trait`.
If `trunc=true` then the distance is truncated to binary value 0 or 1.
196
"""
197
198
199
200
201
202
function get_pairwise_average_isolation(world;trait=1,trunc=false)
        xhist,ttot = get_xhist_mat(world)
        if trunc
                V = truncvar.([xhist[:,i] for i in 1:size(xhist,2)])
        else
                V = var(xhist,dims=1,corrected=false)[:]
203
        end
204
        return sum(V .* (ttot[2:end] .- ttot[1:end-1]))
205
206
207
end

"""
208
    function get_local_pairwise_average_isolation(world,dist,trait=1)
209
210
211
Similar to `get_pairwise_average_isolation`, but the pairwise distance is calculated within demes.
An average of this metrics by deme is return.
"""
212
function get_local_pairwise_average_isolation(world;trait=1,trunc=false)
213
214
        f(a) = get_x(a,1)
        groups = groupby(f,world)
215
        d = get_pairwise_average_isolation.(collect(values(groups)),trait=trait,trunc=trunc)
216
217
        mean(d)
end