Commit eb64136a authored by Victor's avatar Victor
Browse files

debugging hamming metrics

parent 9996f202
...@@ -121,18 +121,18 @@ function get_beta_div(world::World,trait=1) ...@@ -121,18 +121,18 @@ function get_beta_div(world::World,trait=1)
end end
""" """
function get_xhist_mat(world,trait=1,time = 0) $(SIGNATURES)
returns `xhist,ttot`, where `xhist` is a matrix with dimension `lenght(world)` x `length(thist)+1`, returns `xhist,ttot`, where `xhist` is a matrix with dimension `lenght(world)` x `length(thist)+1`,
which consists of geographical history of ancestors at every time step. 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, 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 then the last column of xhist corresponds to actual position of agents
""" """
function get_xhist_mat(world,trait=1,time = 0) function get_xhist_mat(agentarray::Vector{A},trait=1,time = 0) where {A<:AbstractAgent}
thist = vcat(get_thist.(world)...) thist = vcat(get_thist.(agentarray)...)
ttot = sort!(unique(thist)) ttot = sort!(unique(thist))
xhist = zeros(length(world),length(ttot)) xhist = zeros(length(agentarray),length(ttot))
for (i,a) in enumerate(world) for (i,a) in enumerate(agentarray)
_thist = get_thist(a) _thist = get_thist(a)
# Here we check that there is no redundancy of thist because of casting errors # 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 # In the future, we should remove this check, as the type of time has been set to Float64
...@@ -194,8 +194,8 @@ end ...@@ -194,8 +194,8 @@ end
Returns the integrated pairwise squared distance between all agents of `world` wrt `trait`. 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. If `trunc=true` then the distance is truncated to binary value 0 or 1.
""" """
function get_pairwise_average_isolation(world;trait=1,trunc=false) function get_pairwise_average_isolation(world::World;trait=1,trunc=false)
xhist,ttot = get_xhist_mat(world) xhist,ttot = get_xhist_mat(agents(world))
if trunc if trunc
V = truncvar.([xhist[:,i] for i in 1:size(xhist,2)]) V = truncvar.([xhist[:,i] for i in 1:size(xhist,2)])
else else
...@@ -209,9 +209,14 @@ end ...@@ -209,9 +209,14 @@ end
Similar to `get_pairwise_average_isolation`, but the pairwise distance is calculated within demes. Similar to `get_pairwise_average_isolation`, but the pairwise distance is calculated within demes.
An average of this metrics by deme is return. An average of this metrics by deme is return.
""" """
function get_local_pairwise_average_isolation(world;trait=1,trunc=false) function get_local_pairwise_average_isolation(world::World;trait=1,trunc=false)
f(a) = get_x(a,1) f(a) = a[1]
groups = groupby(f,world) groups = groupby(f,agents(world))
d = get_pairwise_average_isolation.(collect(values(groups)),trait=trait,trunc=trunc) 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)
mean(d) mean(d)
end end
cd(@__DIR__)
using Random
Random.seed!(0)
using LightGraphs
using Test
using Revise,ABMEv
using UnPack,JLD2
myspace = (RealSpace{1,Float64}(),)
sigma_K = .9;
sigma_a = .7;
K0 = 1000;
b(X) = gaussian(X[1],0.,sigma_K)
d(X,Y) = gaussian(X[1],Y[1],sigma_a)/K0
D = (1e-2,)
mu = [.1]
NMax = 10000
tend = 1.5
p = Dict{String,Any}();@pack! p = d,b,D,mu,NMax
myagents = [Agent(myspace,(0,),ancestors=true,rates=true) for i in 1:K0]
w0 = World(myagents,myspace,p,0.)
@info "Running simulation with Gillepsie algorithm"
@time sim = run!(w0,Gillepsie(),tend)
@testset "Hamming distances" begin
@test typeof(get_xhist_mat(agents(w0))[1] )<: Array
@test get_pairwise_average_isolation(w0) >0
@test get_local_pairwise_average_isolation(w0) > 0
end
...@@ -4,6 +4,7 @@ using ABMEv, Test, JLD2,Random ...@@ -4,6 +4,7 @@ using ABMEv, Test, JLD2,Random
include("gillepsie.jl") include("gillepsie.jl")
# include("wrightfisher.jl") # include("wrightfisher.jl")
include("metrics.jl") include("metrics.jl")
include("metrics_hamming.jl")
include("utilstest.jl") include("utilstest.jl")
include("simulation.jl") include("simulation.jl")
include("space_agent.jl") include("space_agent.jl")
......
Supports Markdown
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