diff --git a/Project.toml b/Project.toml index 194811487..d31fde164 100644 --- a/Project.toml +++ b/Project.toml @@ -2,7 +2,7 @@ name = "EcologicalNetworks" uuid = "f03a62fe-f8ab-5b77-a061-bb599b765229" authors = ["Timothée Poisot "] repo = "https://github.com/EcoJulia/EcologicalNetworks.jl" -version = "0.5.1" +version = "0.5.2" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" diff --git a/docs/make.jl b/docs/make.jl index 035dfe2c3..0906e2b85 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -30,7 +30,8 @@ makedocs( "Nestedness" => "properties/nestedness.md", "Motifs" => "properties/motifs.md", "Centrality and paths" => "properties/paths.md", - "Overlap and similarity" => "properties/overlap.md" + "Overlap and similarity" => "properties/overlap.md", + "Food web measures" => "properties/foodwebs.md" ], "Advanced information" => [ "Beta-diversity" => "properties/betadiversity.md", diff --git a/docs/src/properties/foodwebs.md b/docs/src/properties/foodwebs.md new file mode 100644 index 000000000..6bcf9b5d0 --- /dev/null +++ b/docs/src/properties/foodwebs.md @@ -0,0 +1,7 @@ +# Food web measures + +```@docs +distance_to_producer +trophic_level +omnivory +``` diff --git a/src/EcologicalNetworks.jl b/src/EcologicalNetworks.jl index f90e928d4..a25ea4e97 100644 --- a/src/EcologicalNetworks.jl +++ b/src/EcologicalNetworks.jl @@ -162,10 +162,7 @@ export KGL01, KGL02, KGL03, KGL04, KGL05, KGL06, KGL07, KGL08, KGL09, KGL10, # Food webs include(joinpath(".", "foodwebs/trophiclevels.jl")) -export fractional_trophic_level, trophic_level - -include(joinpath(".", "foodwebs/omnivory.jl")) -export omnivory +export distance_to_producer, trophic_level, omnivory include(joinpath(".", "information/entropy.jl")) export entropy, make_joint_distribution, mutual_information, conditional_entropy, diff --git a/src/foodwebs/omnivory.jl b/src/foodwebs/omnivory.jl deleted file mode 100644 index 9f6ee28e1..000000000 --- a/src/foodwebs/omnivory.jl +++ /dev/null @@ -1,24 +0,0 @@ -function omnivory(N::T) where {T <: UnipartiteNetwork} - OI = Dict([s => 0.0 for s in species(N)]) - - TL = fractional_trophic_level(N) - k = degree(N; dims=1) - - for sp_i in species(N) - - # Species with no interaction have an omnivory index of 0 - k[sp_i] > 0 || continue - - # For every species, we set its initial omnivory to 0 - oi = 0.0 - for (j, sp_j) in enumerate(species(N)) - # Then for every species it consumes, we ha - tl_diff = (TL[sp_j] - (TL[sp_i]-1.0)).^2.0 - corr = N[sp_i,sp_j]/k[sp_i] - oi += tl_diff * corr - end - OI[sp_i] = oi - end - - return OI -end diff --git a/src/foodwebs/trophiclevels.jl b/src/foodwebs/trophiclevels.jl index b3b23bbb1..3a7ecb739 100644 --- a/src/foodwebs/trophiclevels.jl +++ b/src/foodwebs/trophiclevels.jl @@ -1,24 +1,107 @@ -function fractional_trophic_level(N::T) where {T<:UnipartiteNetwork} - Y = nodiagonal(N) - producers = keys(filter(spedeg -> spedeg.second == 0, degree(Y; dims=1))) - sp = shortest_path(Y) - prod_id = findall(isequal(0), vec(sum(sp; dims=2))) - return Dict(zip(species(Y; dims=1), maximum(sp[:,prod_id]; dims=2).+1)) +""" + distance_to_producer(N::UnipartiteNetwork; f=maximum) + +Returns the distance from a primary producert of every species in the food +web. Primary producers have a trophic level of 1, after that, it's complicated +(see *e.g.* Williams & Martinez 2004, Thompson et al. 2007). + +Post (2002) notes that the trophic level can be obtained from the maximum, +mean, or minimum distance to a producer. Given that consumers may be connected +to more than one producer, one might argue that the *mode* of this distribution +is also relevant. + +In favor of the minimum, one can argue that most energy transfer should happen +along short chains; but imagining a consumer atop a chain of length 5, also +connected directly to a producer, the minimum would give it a trophic level +of 2, hiding its position at the top of the food web. + +In favor of the maximum, one can argue that the higher chains give a better +idea of how far energy coming from the bottom of the food web can go. This +is a strong indication of how *vertically* diverse it is. + +In this package, the keyword `f` defaults to `maximum`. Using any other +function, as long as it accepts an array (of chain distances) and return a +scalar, will work; `minimum` and `Statistics.mean` are natural alternatives, +as is `StatsBase.mode` or `Statistics.median`. + +The chain length from any species to any producer is taken as the shortest +possible path between these two species, as identified by the Dijkstra +algorithm. +""" +function distance_to_producer(N::T; f=maximum) where {T<:UnipartiteNetwork} + Y = nodiagonal(N) + paths = dijkstra(Y) + consumers = collect(keys(filter(p -> iszero(p.second), degree(Y; dims=1)))) + filter!(int -> int.to ∈ consumers, paths) + tl = Dict{eltype(species(Y)),Float64}() + for sp in species(Y) + sp_paths = filter(int -> isequal(sp)(int.from), paths) + chain_lengths = [int.weight for int in sp_paths] + tl[sp] = isempty(chain_lengths) ? 1.0 : f(chain_lengths) + 1.0 + end + return tl end -function trophic_level(N::T) where {T<:UnipartiteNetwork} - TL = fractional_trophic_level(N) - Y = nodiagonal(N) - D = zeros(Float64, size(Y)) - ko = degree_out(Y) +""" + trophic_level(N::UnipartiteNetwork; kwargs...) + +Returns the *fractional* trophic level (after Odum & Heald 1975) +of species in a binary unipartite network. The trophic level is +calculated after Pauly & Christensen (1995), specifically as TLᵢ = 1 + +∑ⱼ(TLⱼ×DCᵢⱼ)/∑ⱼ(DCᵢⱼ), wherein TLᵢ is the trophic level +of species i, and DCᵢⱼ is the proportion of species j in the diet of +species i. Note that this function is calculated on a network where the +diagonal (i.e. self loops) are removed. + +This function uses a *pseudo*-inverse to solve the linear system +described above *iff* the determinant of the diet matrix is 0 (it is +non-invertible). Otherwise, it will use an exact inverse. +""" +function trophic_level(N::UnipartiteNetwork) + # Get the degree as a vector, ordered the same way as species + kout = EcologicalNetworks.degree_out(N) + 𝐤 = [kout[s] for s in species(N)] + # Adjacency matrix to solve the TL + 𝐀 = adjacency(N) + # Diet matrix + 𝐃 = .-(𝐀 ./ 𝐤) + replace!(𝐃, NaN => -0.0) + 𝐃[diagind(𝐃)] .= 1.0 .- 𝐃[diagind(𝐃)] + # Solve with the inverse matrix + 𝐛 = ones(eltype(𝐃), richness(N)) + invfunc = iszero(det(𝐃)) ? pinv : inv + return Dict(zip(species(N), invfunc(𝐃) * 𝐛)) +end + +""" + omnivory(N::T) where {T <: UnipartiteNetwork} + +Returns a vector of length richness(N) with an index of consumption across +trophic levels. Note we explicitly assume that diet contribution is spread +evenly across prey species, in proportion to their degree. Omnivory is +measured based on the output of `trophic_level`. + +#### References +- Christensen, Villy, and Daniel Pauly. "ECOPATH II—a software for balancing + steady-state ecosystem models and calculating network characteristics." + Ecological modelling 61, no. 3-4 (1992): 169-185. + +""" +function omnivory(N::T) where {T<:UnipartiteNetwork} + # Get the trophic level as an array + TL = trophic_level(N) + + # Prepare the omnivory index dictionary + OI = Dict(zip(species(N), fill(0.0, richness(N)))) - # inner loop to avoid dealing with primary producers - for i in 1:richness(Y) - if ko[species(Y)[i]] > 0.0 - D[i,:] = Y[i,:]./ko[species(Y)[i]] + # Loop + for sp in species(N) + if !isempty(N[sp,:]) + TLj = [TL[prey] for prey in N[sp,:]] + # Note that we don't need the normalisation here, because ∑Dᵢⱼ=1 + OI[sp] = sum((TLj .- mean(TLj)).^2.0) + end end - end - # return - return Dict(zip(species(N), 1 .+ D * [TL[s] for s in species(Y)])) + return OI end diff --git a/test/community/foodwebs.jl b/test/community/foodwebs.jl index f1f57f310..836e0783d 100644 --- a/test/community/foodwebs.jl +++ b/test/community/foodwebs.jl @@ -1,27 +1,26 @@ module TestFoodWebs - using Test - using EcologicalNetworks +using Test +using EcologicalNetworks - # test trophic levels on motifs - s1 = unipartitemotifs()[:S1] - s1_ftl = fractional_trophic_level(s1) - s1_tl = trophic_level(s1) - @test maximum(values(s1_ftl)) == 3 - @test maximum(values(s1_tl)) == 3.0 - @test s1_ftl[:a] == 3 - @test s1_ftl[:b] == 2 - @test s1_ftl[:c] == 1 +# test trophic levels on motifs +s1 = unipartitemotifs()[:S1] +s1_ftl = trophic_level(s1) +s1_tl = distance_to_producer(s1) +@test maximum(values(s1_ftl)) ≈ 3 +@test maximum(values(s1_tl)) ≈ 3.0 +@test s1_ftl[:a] ≈ 3 +@test s1_ftl[:b] ≈ 2 +@test s1_ftl[:c] ≈ 1 - s2 = unipartitemotifs()[:S2] - @test minimum(values(fractional_trophic_level(s1))) == 1.0 - @test minimum(values(trophic_level(s1))) == 1.0 +s2 = unipartitemotifs()[:S2] +@test minimum(values(trophic_level(s1))) ≈ 1.0 +@test minimum(values(distance_to_producer(s1))) == 1.0 - # test on a quantitative bipartite network - @test web_of_life("A_HP_002") |> - x -> convert(BinaryNetwork, x) |> x -> convert(UnipartiteNetwork, x) |> - trophic_level |> values |> minimum |> x -> x == 1.0 - @test web_of_life("A_HP_002") |> - x -> convert(BinaryNetwork, x) |> x -> convert(UnipartiteNetwork, x) |> - trophic_level |> values |> maximum |> x -> x == 2.0 +# Test for omnivory +for s in species(s1) + if iszero(EcologicalNetworks.degree_out(s1)[s]) + @test iszero(omnivory(s1)[s]) + end +end end