ABMEv_metrics.jl 6.96 KB
 Victor Boussange committed Jan 17, 2020 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 `````` Victor committed May 07, 2020 45 ``````# TODO: rename this to gamma diversity `````` Victor Boussange committed Jan 17, 2020 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 """ `````` Victor committed Oct 01, 2020 54 ``````function var(world::World;trait=1) `````` Victor committed May 07, 2020 55 `````` xarray = get_x(world,trait) `````` Victor committed Jun 01, 2020 56 `````` return var(xarray,dims=1,corrected=false) `````` Victor Boussange committed Jan 17, 2020 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 committed Oct 01, 2020 66 67 ``````function covgeo(world::World,trait = 0) xarray = Float64.(get_geo(world)) `````` Victor Boussange committed Jan 17, 2020 68 `````` if trait > 0 `````` Victor committed Oct 02, 2020 69 `````` xstd = get_x(world,trait) `````` Victor Boussange committed Jan 17, 2020 70 71 `````` xarray = hcat(xarray,xstd) end `````` Victor committed Jun 01, 2020 72 `````` return cov(xarray,corrected=false) `````` Victor Boussange committed Jan 17, 2020 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 committed Oct 01, 2020 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)) `````` Victor Boussange committed Jan 17, 2020 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 `````` Victor committed May 07, 2020 91 ``````""" `````` Victor committed May 08, 2020 92 `````` get_alpha_div(world::Array{U,1},t::Number,trait=1) where U <: Union{Missing,Agent} `````` Victor committed Oct 01, 2020 93 ``````Mean of the local variance of `trait` per patch. If trait=0, we get the geotrait `````` Victor committed May 07, 2020 94 95 ``````# Arguments """ `````` Victor committed Oct 01, 2020 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 committed May 07, 2020 99 `````` if trait == 0 `````` Victor committed Oct 01, 2020 100 `````` return mean([var(Float64.([x[end] for x in xp]),corrected=false) for xp in values(g)]) `````` Victor committed May 07, 2020 101 `````` else `````` Victor committed Oct 01, 2020 102 `````` return mean([var(Float64.([x[trait] for x in xp]),corrected=false) for xp in values(g)]) `````` Victor committed May 07, 2020 103 `````` end `````` Victor committed May 07, 2020 104 105 106 ``````end """ `````` Victor committed May 08, 2020 107 `````` get_beta_div(world::Array{U,1},t::Number,trait=1) where U <: Union{Missing,Agent} `````` Victor committed May 07, 2020 108 ``````Variance of the mean of `trait` per patch `````` Victor committed May 07, 2020 109 110 ``````# Arguments """ `````` Victor committed Oct 01, 2020 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 committed May 07, 2020 114 `````` if trait == 0 `````` Victor committed May 08, 2020 115 `````` # need to convert to Float64, otherwise infinite variance `````` Victor committed Oct 01, 2020 116 `````` sbar_i = [mean(Float64.([x[end] for x in xp])) for xp in values(g)] `````` Victor committed May 07, 2020 117 `````` else `````` Victor committed Oct 01, 2020 118 `````` sbar_i = [mean(Float64.([x[trait] for x in xp])) for xp in values(g)] `````` Victor committed May 07, 2020 119 `````` end `````` Victor committed Jun 01, 2020 120 `````` return var(sbar_i,corrected=false) `````` Victor committed May 07, 2020 121 ``````end `````` Victor committed Sep 11, 2020 122 123 `````` """ `````` Victor committed Sep 20, 2020 124 `````` function get_xhist_mat(world,trait=1,time = 0) `````` Victor committed Sep 20, 2020 125 ``````returns `xhist,ttot`, where `xhist` is a matrix with dimension `lenght(world)` x `length(thist)+1`, `````` Victor committed Sep 20, 2020 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 `````` Victor committed Sep 11, 2020 129 ``````""" `````` Victor committed Sep 20, 2020 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) `````` Victor committed Sep 20, 2020 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 `````` Victor committed Sep 20, 2020 146 `````` _xhist = get_xhist(a,trait) `````` Victor committed Sep 20, 2020 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] `````` Victor committed Sep 20, 2020 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) `````` Victor committed Sep 20, 2020 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])) `````` Victor committed Sep 11, 2020 179 180 ``````end `````` Victor committed Sep 21, 2020 181 182 183 184 ``````""" function truncvar(X::AbstractArray) Returns the truncated variance of `X`. """ `````` Victor committed Sep 20, 2020 185 `````` `````` Victor committed Sep 21, 2020 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 `````` Victor committed Sep 20, 2020 191 `````` `````` Victor committed Sep 11, 2020 192 ``````""" `````` Victor committed Sep 21, 2020 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. `````` Victor committed Sep 11, 2020 196 ``````""" `````` Victor committed Sep 21, 2020 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)[:] `````` Victor committed Sep 11, 2020 203 `````` end `````` Victor committed Sep 21, 2020 204 `````` return sum(V .* (ttot[2:end] .- ttot[1:end-1])) `````` Victor committed Sep 11, 2020 205 206 207 ``````end """ `````` Victor committed Sep 20, 2020 208 `````` function get_local_pairwise_average_isolation(world,dist,trait=1) `````` Victor committed Sep 11, 2020 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. """ `````` Victor committed Sep 21, 2020 212 ``````function get_local_pairwise_average_isolation(world;trait=1,trunc=false) `````` Victor committed Sep 11, 2020 213 214 `````` f(a) = get_x(a,1) groups = groupby(f,world) `````` Victor committed Sep 21, 2020 215 `````` d = get_pairwise_average_isolation.(collect(values(groups)),trait=trait,trunc=trunc) `````` Victor committed Sep 11, 2020 216 217 `````` mean(d) end``````