Merge pull request #185 from EcoJulia/master

RDPG

docs/src/random/null.md

@@ -24,11 +24,23 @@ null2
 null3
 ```
 
+### Random Dot Product Graph model
+
+The interaction probability in a Random Dot Product Graphs model are given by the dot product of species representations in two metric spaces of a given dimension (one describing species as consumers, one as producers -- or predators and preys).
+
+```@docs
+RDPG
+```
+
+The Random Dot Product Graph spaces are computed via a truncated Singular Value Decomposition of the food web adjacency matrix.
+
+```@docs
+svd_truncated
+```
+
 ## Shuffle interactions
 
 ```@docs
 shuffle!
 shuffle
 ```
-
-

src/EcologicalNetworks.jl

@@ -75,6 +75,10 @@ export isdegenerate, simplify, simplify!
    #, species_has_no_successors, species_has_no_predecessors,
    #species_is_free, free_species
 
+# SVD
+include(joinpath(".", "svd", "svd.jl"))
+export rdpg
+
 # Random networks and permutations
 include(joinpath(".", "rand/draws.jl"))
 include(joinpath(".", "rand/sample.jl"))
@@ -84,6 +88,8 @@ include(joinpath(".", "rand/null.jl"))
 include(joinpath(".", "rand/linfilter.jl"))
 export linearfilter, linearfilterzoo
 export null1, null2, null3, null4
+include(joinpath(".", "rand/RDPG.jl"))
+export RDPG
 
 # Random networks from structural models
 include(joinpath(".", "structuralmodels/cascademodel.jl"))

src/rand/RDPG.jl

@@ -0,0 +1,28 @@
+"""
+    RDPG(N::BinaryNetwork; rank::Integer=3)
+
+Given a binary network `N`, `RDPG(N)` returns a probabilistic network with the
+same number of species, where every interaction happens with a probability equal
+to the dot product of species representation in the network `N`'s RDPG space of
+rank `rank`.
+
+Because the pairwise dot product obtained by the matrix multiplication of the
+two spaces `Left * Right` are not granted to be bounded between 0 and 1 (for
+numerical and theoric reasons), we bound the entries to be in the `[0,1]` range.
+
+
+#### References
+
+Dalla Riva, G.V. and Stouffer, D.B., 2016. Exploring the evolutionary signature
+of food webs' backbones using functional traits. Oikos, 125(4), pp.446-456.
+https://doi.org/10.1111/oik.02305
+
+"""
+function RDPG(N::T; rank::Integer=3) where {T<:BinaryNetwork}
+  L, R = svd_truncated(N; rank=rank)
+  ardpg = L * R
+  ardp[ardpg .< 0.] .= 0.
+  ardpg[ardpg .> 1.] .= 1.
+  ReturnType = typeof(N) <: AbstractBipartiteNetwork ? BipartiteProbabilisticNetwork : UnipartiteProbabilisticNetwork
+  return ReturnType(ardpg, EcologicalNetworks.species_objects(N)...)
+end

src/svd/svd.jl

@@ -0,0 +1,45 @@
+import LinearAlgebra: svd, mul!, diagm
+
+"""
+    LinearAlgebra.svd(N::T) where {T <: AbstractEcologicalNetwork}
+
+SVD of a network (i.e. SVD of the adjacency matrix)
+"""
+LinearAlgebra.svd(N::T) where {T <: AbstractEcologicalNetwork} = svd(adjacency(N))
+
+"""
+    svd_truncated(N::BinaryNetwork, rnk::Integer=3)
+
+Given a binary network `N` which adjacency matrix `A` is of dimension `n × m`,
+`svd_truncated(A)` returns two matrices, `Left` and `Right`, with dimensions
+respectively `n × rank` and `rank × m`, corresponding to the species
+representation in the network `N`'s RDPG space of rank `rank`.
+
+The singular value decomposition is computed using `LinearAlgebra`'s `svd`,
+obtaining
+
+`A = U * Diagonal(S) * V = U * Diagonal(√S) * Diagonal(√S) * V`.
+
+The truncation preserves the first `rank` columns of `U * Diagonal(√S)` and the
+first `rank` rows `Diagonal(√S) * V`.
+
+We have that, `A ≃ Left * Right` (and the approximation is optimal in a
+specified meaning).
+
+#### References
+
+Dalla Riva, G.V. and Stouffer, D.B., 2016. Exploring the evolutionary signature
+of food webs' backbones using functional traits. Oikos, 125(4), pp.446-456.
+https://doi.org/10.1111/oik.02305
+
+"""
+function rdpg(N::T, rnk::Integer=3) where {T<:BinaryNetwork}
+    U, singular_values, V = svd(N)
+    left_space = similar(U)
+    right_space = similar(V')
+    mul!(left_space, U, diagm(sqrt.(singular_values)))
+    mul!(right_space, diagm(sqrt.(singular_values)), V')
+    left_space = left_space[:,1:rnk]
+    right_space = right_space[1:rnk,:]
+    return left_space, right_space
+end