Commit 4f836c2e authored by Victor's avatar Victor
Browse files

added get_local_abundance and modified get_alpha_div

parent 1f6a4ad1
Pipeline #80583 failed with stage
in 31 minutes and 13 seconds
......@@ -34,7 +34,7 @@ module ABMEv
updateDeathEvent!#,runWorld_G!,runWorld_WF!,
export Simulation,add_entry!,get_tend,get_size,get_tspan,get_world,get_xnt
export H_discrete,findclusters,var,covgeo,hamming,get_beta_div, get_alpha_div,
get_dist_hist,get_pairwise_average_isolation,
get_local_abundance,get_dist_hist,get_pairwise_average_isolation,
get_local_pairwise_average_isolation,
truncvar,get_xhist_mat
export update_afterbirth_std!,update_afterdeath_std!
......
......@@ -56,6 +56,7 @@ function var(world::World;trait=1)
xarray = get_x(world,trait)
if trait > 0
if typeof(space(world)[trait]) <: GraphSpace
# not working
fiedlervec = eigs(laplacian_matrix(space(world)[trait].g),nev=2,which=:SM)[2][:,2]
return mean(fiedlervec[xarray].^2) - mean(fiedlervec[xarray])^2
end
......@@ -116,13 +117,37 @@ end
get_alpha_div(world::Array{U,1},t::Number,trait=1) where U <: Union{Missing,Agent}
Mean of the local variance of `trait` per patch.
If trait=0, we get the mean of the local variance of the geotrait
If average = false, returns the alpha div for each patch, ordered by vertices
"""
function get_alpha_div(world::World,trait=1)
function get_alpha_div(world::World,trait=1,average=true)
g = groupby(a->a[1],agents(world))
# here the second mean is here when subspace is multidimensional
v = [var(World(subw,space(world),parameters(world)),trait=trait) for subw in values(g)]
# we sort by index of vertices
v = [var(World(g[i],space(world),parameters(world)),trait=trait) for i in sort(collect(keys(g)))]
h = vcat(v...)
return mean(h)
if average
return mean(h)
else
return mean(h,dims=2)
end
end
"""
get_alpha_div(world::Array{U,1},t::Number,trait=1) where U <: Union{Missing,Agent}
Mean of the local variance of `trait` per patch.
If trait=0, we get the mean of the local variance of the geotrait
If average = false, returns the alpha div for each patch, ordered by vertices
"""
function get_local_abundance(world::World,average=true)
g = groupby(a->a[1],agents(world))
# here the second mean is here when subspace is multidimensional
# we sort by index of vertices
a = [length(g[i]) for i in sort(collect(keys(g)))]
if average
return mean(a)
else
return a
end
end
"""
......
......@@ -4,14 +4,14 @@ using Revise,ABMEv
using UnPack
myspace1 = (RealSpace{1,Float64}(),)
myspace2 = (RealSpace{1,Float64}(),RealSpace{1,Float64}())
myspace3 = (DiscreteSegment(Int16(1),Int16(10)),RealSpace{1,Float64}())
myspace3 = (DiscreteSegment(Int16(1),Int16(10)),RealSpace{2,Float64}())
g = SimpleGraph(1000,4000)
myspace4 = (RealSpace{1,Float64}(),GraphSpace(g),)
K0 = 1000; σ = 1e-1
a1 = [Agent(myspace1, (σ,) .* randn() .- .5,ancestors=true) for i in 1:K0]
a2 = [Agent(myspace2,tuple((σ, σ) .* randn(2) .- .5...),ancestors=true) for i in 1:K0]
a3 = [Agent(myspace3, (rand(Int16.(1:10)), 1e-1* randn() + 5.5 ),ancestors=true) for i in 1:K0]
a3 = [Agent(myspace3, (rand(Int16.(1:10)), tuple((1e-1.* randn(2) .+ 5.5)... )),ancestors=true) for i in 1:K0]
a4 = [Agent(myspace4, (1.,rand(Int64(1):Int64(1000)),),ancestors=true) for i in 1:K0]
D = (1.,);
......@@ -20,7 +20,7 @@ NMax = 1000
p1 = Dict{String,Any}();@pack! p1 = D,mu,NMax
D = (1.,1.);
p2 = Dict{String,Any}();@pack! p2 = D,mu,NMax
D = (Int16(0.),0.)
D = (Int16(0.),(0.,0.))
p3 = Dict{String,Any}();@pack! p3 = D,mu,NMax
D = (0.,0.)
p4 = Dict{String,Any}();@pack! p4 = D,mu,NMax
......@@ -36,9 +36,7 @@ w4 = World(a4,myspace4,p4)
@test first(var(w1)) (σ).^2 atol=0.001
@test first(var(w2,trait=2)) (σ).^2 atol=0.001
@test first(var(w4,trait=1)) < Inf
@test first(var(w4,trait=2)) < Inf
# @test first(var(w4,trait=2)) < Inf
end
## testing covgeo
......@@ -60,7 +58,14 @@ w4 = World(a4,myspace4,p4)
@test abs(α) < Inf
α = get_alpha_div(w4,1);
@test abs(α) < Inf
α = get_alpha_div(w3,2,false)
@test length(α) == 10
end
@testset "Abundance" begin
@test isapprox(get_local_abundance(w3), K0/10,atol = 10)
N = get_local_abundance(w3,false)
@test length(N) == 10
end
@testset "Beta diversity" begin
β = get_beta_div(w3,2);
@test abs(β) < Inf
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment