added hamming distance for x history metrics, and changed structure of agents...

added hamming distance for x history metrics, and changed structure of agents so that t_hist is now a Float64 (otherwise, conflicts happen because dt can be less that Float16 error)

Project.toml

 name = "ABMEv"
 uuid = "837ac870-fb52-4b0c-9a0e-030f2f36f5ed"
 authors = ["Victor Boussange "]
-version = "0.2.1"
+version = "0.2.2"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"

src/ABMEv.jl

@@ -16,11 +16,13 @@ module ABMEv
     @reexport using Distributions, DataFrames
     export update_rates!
     export MixedAgent,StdAgent,Agent,get_fitness,get_x,get_dim,get_nancestors,get_xarray,get_xhist,
-        get_geo,get_b,get_d,increment_x!,get_inc_reflected,world2df,
+        get_thist,get_geo,get_b,get_d,increment_x!,get_inc_reflected,world2df,
         split_move,split_merge_move,tin,new_world_G
     export copy,runWorld_store_WF,runWorld_store_G,clean_world #,runWorld_G!,runWorld_WF!,
-    export H_discrete,findclusters,var,covgeo,hamming,get_beta_div, get_alpha_div
+    export H_discrete,findclusters,var,covgeo,hamming,get_beta_div, get_alpha_div,
+        get_hamming_dist_hist,get_pairwise_average_isolation,
+        get_local_pairwise_average_isolation,
     export update_afterbirth_std!,update_afterdeath_std!
     export generalised_gaussian,gaussian,ma,geomsmooth,arithsmooth,eth_grad_std,
-        DiversityFunction,geomsmooth2D,arithsmooth2D,interpolate_df
+        DiversityFunction,geomsmooth2D,arithsmooth2D,interpolate_df,groupby
 end

src/ABMEv_Agent.jl

@@ -5,7 +5,7 @@ mutable struct Agent{T,U}
     # history of traits for geotraits
     x_history::Array{U}
     # birth time of ancestors
-    t_history::Array{U,1}
+    t_history::Array{Float64,1}
     # death rate
     d::Float64
     #birth rate
@@ -49,6 +49,7 @@ get_x(a::Agent,t::Number,i::Integer) = i > 0 ? a.x_history[Int(i),end] : get_geo
 get_x(a::Agent,i::Integer) = a.x_history[Int(i),end]
 get_xhist(a::Agent,i::Number) = a.x_history[Int(i),:]
 get_xhist(a::Agent) = a.x_history
+get_thist(a::Agent) = a.t_history
 get_d(a::Agent) = a.d
 get_b(a::Agent) = a.b
 get_fitness(a::Agent) = a.b - a.d

src/ABMEv_metrics.jl

@@ -122,3 +122,66 @@ function get_beta_div(world::Array{U,1},t::Number,trait=1) where U <: Union{Miss
     end
     return var(sbar_i,corrected=false)
 end
+
+"""
+    function get_hamming_dist_hist(a1,a2,time,trait=1)
+Returns a metric corresponding to the exact time agents have been isolated to each other
+"""
+function get_hamming_dist_hist(a1,a2,trait=1,time = 0)
+        thist = vcat(get_thist(a1),get_thist(a2))
+        ttot = sort!(unique(thist))
+        xhist = zeros(2,length(ttot))
+        for (i,a) in enumerate([a1,a2])
+                _thist = get_thist(a)
+                _xhist = get_xhist(a,trait)
+                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]
+        end
+        tt = 0
+        time > ttot[end] ? ttot = vcat(ttot,time) : tt = 1
+        return sum((xhist[1,1:end-tt] .!= xhist[2,1:end-tt])[:] .* (ttot[2:end] .- ttot[1:end-1]))
+        # return xhist
+        # return (thist[2:end] .- thist[1:end-1])
+end
+
+"""
+    function get_pairwise_average_isolation(world,trait=1,time = 0)
+Returns the average pairwise isolation by time between all agents of `world`,
+using metrics `get_hamming_dist_hist`
+"""
+function get_pairwise_average_isolation(world,trait=1,time = 0)
+        N = size(world,1)
+        D = zeros(N)
+        # Here you should do a shared array to compute in parallel
+        Threads.@threads for i in 1:(N-1)
+                for j in i+1:N
+                            C = get_hamming_dist_hist(world[i],world[j],trait,time)
+                            D[i] += C
+                            D[j] += C
+                end
+        end
+        return sum(D) / N^2
+end
+
+"""
+    function get_local_pairwise_average_isolation(world,trait=1,time = 0)
+Similar to `get_pairwise_average_isolation`, but the pairwise distance is calculated within demes.
+An average of this metrics by deme is return.
+"""
+function get_local_pairwise_average_isolation(world,trait=1,time = 0)
+        f(a) = get_x(a,1)
+        groups = groupby(f,world)
+        d = get_pairwise_average_isolation.(collect(values(groups)),trait,time)
+        mean(d)
+end

src/ABMEv_utils.jl

@@ -46,7 +46,7 @@ const eth_grad_small = cgrad([colorant"#1F407A", RGB(0.671,0.851,0.914),RGB(1.0,
 const eth_grad_std = cgrad([colorant"#1F407A", RGB(0.671,0.851,0.914),RGB(1.0,1.0,0.749), RGB(0.992,0.682,0.38),RGB(0.647,0.0,0.149),],[.0,1.])
 
 # This is all about interpolations
-using Interpolations
+import Interpolations:interpolate,Gridded,Linear
 struct DiversityFunction
     x
     y
@@ -94,3 +94,16 @@ function interpolate_df(df,xlab,ylab,zlab)
     A = reshape(df[zlab],length(xa),length(ya))
     return DiversityFunction(xa,ya,interpolate((xa,ya),A,Gridded(Linear())))
 end
+
+import DataFrames.groupby
+"""
+    function groupby(f, list::Array)
+returns a dictionary that group `list` elements by value of function `f`
+"""
+groupby(f, list::Array) = begin
+  groups = Dict()
+  for v in list
+    push!(get!(groups, f(v), []), v)
+  end
+  groups
+end

test/metrics.jl

@@ -33,7 +33,11 @@ p = Dict("mu" => [1. 1.],"D" => [0. 0.], "nodes" =>10 )
            β = get_beta_div(agentsd,1.0,2);
            @test abs(β) < Inf
        end
-
+       @testset "Isolation by history hamming distance" begin
+           a1 = Agent{StdAgent,Float64}([1 2 3],[0,1,4],0.,1.)
+           a2 = Agent{StdAgent,Float64}([1 10 3 10],[0,3,4,5],0.,1.)
+           @test get_hamming_dist_hist(a1,a2,1) ≈ 3.0
+       end
 end
 
 @testset "Geotrait computation" begin