Updating metrics of diversity to take into account graph structure. This will...

Updating metrics of diversity to take into account graph structure. This will have to change in the future.

Manifest.toml

@@ -18,6 +18,18 @@ git-tree-sha1 = "2b6845cea546604fb4dca4e31414a6a59d39ddcd"
 uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
 version = "0.0.4"
 
+[[Arpack]]
+deps = ["Arpack_jll", "Libdl", "LinearAlgebra"]
+git-tree-sha1 = "2ff92b71ba1747c5fdd541f8fc87736d82f40ec9"
+uuid = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+version = "0.4.0"
+
+[[Arpack_jll]]
+deps = ["Libdl", "OpenBLAS_jll", "Pkg"]
+git-tree-sha1 = "e214a9b9bd1b4e1b4f15b22c0994862b66af7ff7"
+uuid = "68821587-b530-5797-8361-c406ea357684"
+version = "3.5.0+3"
+
 [[Artifacts]]
 deps = ["Pkg"]
 git-tree-sha1 = "c30985d8821e0cd73870b17b0ed0ce6dc44cb744"
@@ -170,9 +182,9 @@ version = "0.4.0"
 
 [[FFMPEG_jll]]
 deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "LibVPX_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
-git-tree-sha1 = "ef1fb99ef8f4727dd9ea46fc4c10920a955a8162"
+git-tree-sha1 = "13a934b9e74a8722bf1786c989de346a9602e695"
 uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
-version = "4.3.1+3"
+version = "4.3.1+2"
 
 [[FFTW]]
 deps = ["AbstractFFTs", "FFTW_jll", "IntelOpenMP_jll", "Libdl", "LinearAlgebra", "MKL_jll", "Reexport"]
@@ -432,6 +444,12 @@ git-tree-sha1 = "a42c0f138b9ebe8b58eba2271c5053773bde52d0"
 uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
 version = "1.3.4+2"
 
+[[OpenBLAS_jll]]
+deps = ["CompilerSupportLibraries_jll", "Libdl", "Pkg"]
+git-tree-sha1 = "5fae4d1510bdcf7768cc951878b8aa48666c58a8"
+uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
+version = "0.3.10+0"
+
 [[OpenSSL_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
 git-tree-sha1 = "71bbbc616a1d710879f5a1021bcba65ffba6ce58"

Project.toml

@@ -4,6 +4,7 @@ authors = ["Victor Boussange "]
 version = "1.1.0"
 
 [deps]
+Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"

README.md

 # ABMEv.jl
 [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://vboussange.github.io/ABMEv.jl/stable)
 [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://vboussange.github.io/ABMEv.jl/dev)
+[![pipeline status](https://gitlab.ethz.ch/bvictor/abmev/badges/master/pipeline.svg)](https://gitlab.ethz.ch/bvictor/abmev/-/commits/master)
+[![coverage report](https://gitlab.ethz.ch/bvictor/abmev/badges/master/coverage.svg)](https://gitlab.ethz.ch/bvictor/abmev/-/commits/master)
 
 <div align="center"> <img
 src="docs/src/assets/abmev_1d.png"

src/ABMEv.jl

@@ -5,6 +5,7 @@ module ABMEv
     using LightGraphs
     using UnPack
     using DocStringExtensions
+    using Arpack
 
     include("ABMEv_Space.jl")
     include("ABMEv_Agent.jl")

src/ABMEv_metrics.jl

@@ -44,33 +44,46 @@ end
 import Statistics:var,mean
 # TODO: rename this to gamma diversity
 """
-    var(world::Array{Agent};trait=:)
+    function var(world::World;trait=1)
+Return the variance of the `world`'s `trait` distribution.
 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
-
+# Notes
+For now, the variance of a `trait` defined on a `GraphSpace` is calculated
+thanks to the Fiedler vector (cf `https://mathworld.wolfram.com/FiedlerVector.html`)
 """
 function var(world::World;trait=1)
-    xarray = Float64.(get_x(world,trait))
-    return var(xarray,dims=1,corrected=false)
+    xarray = get_x(world,trait)
+    if trait > 0 && typeof(space(world)[trait]) <: GraphSpace
+        fiedlervec = eigs(laplacian_matrix(space(world)[trait].g),nev=2,which=:SM)[2][:,2]
+        return mean(fiedlervec[xarray].^2) - mean(fiedlervec[xarray])^2
+    else
+        return var(Float64.(xarray),dims=1,corrected=false)
+    end
 end
 
 """
-function mean(world::World;trait=1)
-
+    function mean(world::World;trait=1)
+Returns the mean of the `world`'s `trait` distribution.
+If trait = 0, returns the variance of the geotrait,
 """
 function mean(world::World;trait=1)
-    xarray = Float64.(get_x(world,trait))
-    return mean(xarray,dims=1)
+    xarray = get_x(world,trait)
+    if trait > 0 && typeof(space(world)[trait]) <: GraphSpace
+        fiedlervec = eigs(laplacian_matrix(space(world)[trait].g),nev=2,which=:SM)[2][:,2]
+        return mean(fiedlervec[xarray])
+    else
+        return mean(Float64.(xarray),dims=1)
+    end
 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
-
+# Notes
+This might be deprecated in the future
 """
 function covgeo(world::World,trait = 0)
     xarray = Float64.(get_geo(world))
@@ -99,19 +112,15 @@ function hamming(world::World) where T <: Int
 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 geotrait
-# Arguments
+Mean of the local variance of `trait` per patch.
+If trait=0, we get the mean of the local variance of the geotrait
 """
 function get_alpha_div(world::World,trait=1)
     g = groupby(a->a[1],agents(world))
-    if trait == 0
-        return mean([var(Float64.(get_geo(World(subw,space(world),parameters(world)))),corrected=false) for subw in values(g)])
-    else
-        # 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)]
-        h = vcat(v...)
-        return mean(h)
-    end
+    # 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)]
+    h = vcat(v...)
+    return mean(h)
 end
 
 """
@@ -121,15 +130,9 @@ Variance of the mean of `trait` per patch
 """
 function get_beta_div(world::World,trait=1)
     g = groupby(a->a[1],agents(world))
-    if trait == 0
-        # need to convert to Float64, otherwise infinite variance
-        sbar_i = [mean(Float64.(get_geo(World(subw,space(world),parameters(world))))) for subw in values(g)]
-        return var(sbar_i,corrected=false)
-    else
-        m = [mean(World(subw,space(world),parameters(world)),trait=trait) for subw in values(g)]
-        h=vcat(m...)
-        return mean(var(h,dims=1,corrected=false))
-    end
+    m = [mean(World(subw,space(world),parameters(world)),trait=trait) for subw in values(g)]
+    h=vcat(m...)
+    return mean(var(h,dims=1,corrected=false))
 end
 
 """

src/ABMEv_utils.jl

@@ -40,6 +40,16 @@ function geomsmooth(x,smooth)
     return [prod(x[i-smooth + 1:i])^(1/smooth) for i in smooth:length(x)]
 end
 
+"""
+    function geomsmooth(x,y,smooth)
+-`x` is the shifted x-axis vector, due to smoothing
+-`y` is the smoothed value
+Returning a tuple
+"""
+function geomsmooth(x,y,smooth)
+    idx = Int((smooth-1)/2)+1:length(x)- Int((smooth-1)/2)
+    return x[idx],[prod(y[i-smooth + 1:i])^(1/smooth) for i in smooth:length(y)]
+end
 """
     function arithsmooth(x,smooth)
 arithmetic smoothing
@@ -48,9 +58,22 @@ arithmetic smoothing
 function arithsmooth(x,smooth)
     return [sum(x[i-smooth+1:i])/smooth for i in smooth:length(x)]
 end
+"""
+    function arithsmooth(x,y,smooth)
+-`x` is the shifted x-axis vector, due to smoothing
+-`y` is the smoothed value
+Returning a tuple
+"""
+function arithsmooth(x,y,smooth)
+    idx = Int((smooth-1)/2)+1:length(x)- Int((smooth-1)/2)
+    return x[idx],[sum(y[i-smooth+1:i])/smooth for i in smooth:length(y)]
+end
 
 # This is all about interpolations
 import Interpolations:interpolate,Gridded,Linear
+"""
+$(TYPEDEF)
+"""
 struct DiversityFunction
     x
     y
@@ -99,7 +122,7 @@ function interpolate_df(df,xlab,ylab,zlab)
     if length(xa) > 1 && length(ya) > 1
         return DiversityFunction(xa,ya,interpolate((xa,ya),A,Gridded(Linear())))
     else
-        return DiversityFunction(xa,ya,A)
+        return DiversityFunction(xa[:],ya[:],A[:])
     end
 end
 

test/metrics.jl

@@ -5,10 +5,15 @@ using UnPack
 myspace1 = (RealSpace{1,Float64}(),)
 myspace2 = (RealSpace{1,Float64}(),RealSpace{1,Float64}())
 myspace3 = (DiscreteSegment(Int16(1),Int16(10)),RealSpace{1,Float64}())
+g = SimpleGraph(1000,4000)
+myspace4 = (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]
+a4 = [Agent(myspace4, (rand(Int64(1):Int64(1000)),),ancestors=true) for i in 1:K0]
+
 D = (1.,);
 mu = [1.,1.]
 NMax = 1000
@@ -17,16 +22,21 @@ D = (1.,1.);
 p2 = Dict{String,Any}();@pack! p2 = D,mu,NMax
 D = (Int16(0.),0.)
 p3 = Dict{String,Any}();@pack! p3 = D,mu,NMax
+D = (0.)
+p4 = Dict{String,Any}();@pack! p4 = D,mu,NMax
 
 w1 = World(a1,myspace1,p1)
 w2 = World(a2,myspace2,p2)
 w3 = World(a3,myspace3,p3)
+w4 = World(a4,myspace4,p4)
 
 ## testing variance
 @testset "Testing metrics" begin
     @testset "var" begin
         @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
+
     end
 
     ## testing covgeo
@@ -46,10 +56,14 @@ w3 = World(a3,myspace3,p3)
       @testset "Alpha diversity" begin
           α = get_alpha_div(w3,2);
           @test abs(α) < Inf
+          α = get_alpha_div(w4,1);
+          @test abs(α) < Inf
        end
        @testset "Beta diversity" begin
            β = get_beta_div(w3,2);
            @test abs(β) < Inf
+           β = get_beta_div(w4,1);
+           @test abs(β) < Inf
        end
        @testset "Isolation by history - hamming distance" begin
            a1 = Agent((DiscreteSegment(1,10),),[(1,),(2,),(3,)],[0,1,4],ancestors=true)