ABMEv_metrics.jl 7.3 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 = Float64.(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
69
        xstd = get_x(world,trait)
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
function get_alpha_div(world::World,trait=1)
97
    g = groupby(a->a[1],agents(world))
Victor's avatar
Victor committed
98
    if trait == 0
99
        return mean([var(Float64.(get_geo(World(subw,space(world),parameters(world)))),corrected=false) for subw in values(g)])
Victor's avatar
Victor committed
100
    else
101
102
        # here the second mean is here when subspace is multidimensional
        return mean([mean(var(Float64.(get_x(World(subw,space(world),parameters(world)),trait)),corrected=false)) for subw 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
function get_beta_div(world::World,trait=1)
112
    g = groupby(a->a[1],agents(world))
Victor's avatar
Victor committed
113
    if trait == 0
114
        # need to convert to Float64, otherwise infinite variance
115
        sbar_i = [mean(Float64.(get_geo(World(subw,space(world),parameters(world))))) for subw in values(g)]
Victor's avatar
Victor committed
116
    else
117
        sbar_i = [mean(Float64.(get_x(World(subw,space(world),parameters(world)),trait))) for subw in values(g)]
Victor's avatar
Victor committed
118
    end
Victor's avatar
bug fix    
Victor committed
119
    return var(sbar_i,corrected=false)
120
end
121
122

"""
Victor's avatar
Victor committed
123
$(SIGNATURES)
124
returns `xhist,ttot`, where `xhist` is a matrix with dimension `lenght(world)` x `length(thist)+1`,
125
126
127
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
128
"""
129

Victor's avatar
Victor committed
130
131
function get_xhist_mat(agentarray::Vector{A},trait=1,time = 0) where {A<:AbstractAgent}
        thist = vcat(get_thist.(agentarray)...)
132
        ttot = sort!(unique(thist))
Victor's avatar
Victor committed
133
134
        xhist = zeros(length(agentarray),length(ttot))
        for (i,a) in enumerate(agentarray)
135
136
137
138
139
140
141
142
143
144
            _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
145
                _xhist = get_xhist(a,trait)
146
147
148
149
150
151
152
153
154
155
156
157
158
159
            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]
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
        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)
176
177
        xhist,ttot = get_xhist_mat([a1,a2],trait,time)
        return sum(dist.(xhist[1,:],xhist[2,:]) .* (ttot[2:end] .- ttot[1:end-1]))
178
179
end

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

185
186
187
188
189
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
190

191
"""
192
193
194
    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.
195
"""
Victor's avatar
Victor committed
196
197
function get_pairwise_average_isolation(world::World;trait=1,trunc=false)
        xhist,ttot = get_xhist_mat(agents(world))
198
199
200
201
        if trunc
                V = truncvar.([xhist[:,i] for i in 1:size(xhist,2)])
        else
                V = var(xhist,dims=1,corrected=false)[:]
202
        end
203
        return sum(V .* (ttot[2:end] .- ttot[1:end-1]))
204
205
206
end

"""
207
    function get_local_pairwise_average_isolation(world,dist,trait=1)
208
209
210
Similar to `get_pairwise_average_isolation`, but the pairwise distance is calculated within demes.
An average of this metrics by deme is return.
"""
Victor's avatar
Victor committed
211
212
213
214
215
216
217
218
219
function get_local_pairwise_average_isolation(world::World;trait=1,trunc=false)
        f(a) = a[1]
        groups = groupby(f,agents(world))
        smallworlds = []
        for v in values(groups)
            myagents = [a for a in v]
            push!(smallworlds,World(myagents,world.space,world.p,world.t))
        end
        d = get_pairwise_average_isolation.(smallworlds,trait=trait,trunc=trunc)
220
221
        mean(d)
end