Skip to content

Commit

Permalink
🌪️#2: reframe diversity indices and tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@alaindanet
alaindanet committed Mar 8, 2023
1 parent 0e869f7 commit 07e68a9
Show file tree
Hide file tree
Showing 3 changed files with 124 additions and 98 deletions.
127 changes: 53 additions & 74 deletions src/measures/functioning.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,11 +40,7 @@ julia> even = evenness(sol);
"""
function richness(solution; threshold = eps(), kwargs...)
measure_on = extract_last_timesteps(solution; kwargs...)

rich = [
richness(vec(measure_on[:, i]); threshold = threshold) for
i in 1:size(measure_on, 2)
]
rich = richness.(eachcol(measure_on), threshold = threshold)
mean(rich)
end

Expand Down Expand Up @@ -113,12 +109,10 @@ species_persistence(n::AbstractVector; threshold = eps()) =
"""
biomass(solution; kwargs...)
Returns a named tuple of total and species biomass, averaged over the last `last` timesteps.
Returns a named tuple of total and species biomass, averaged over the `last` timesteps.
# Arguments
- `solution`, `last`: See [`richness`](@ref) doc.
kwargs... arguments are dispatched to [`extract_last_timesteps`](@ref). See
[`extract_last_timesteps`](@ref) for the argument details.
Expand Down Expand Up @@ -161,7 +155,8 @@ biomass(vec::AbstractVector;) = (total = mean(vec), sp = mean(vec))
"""
shannon_diversity(solution; threshold, kwargs...)
Computes the Shannon entropy index, i.e. the first Hill number.
Computes the average Shannon entropy index, i.e. the first Hill number,
over the `last` timesteps.
kwargs... arguments are dispatched to [`extract_last_timesteps`](@ref). See
[`extract_last_timesteps`](@ref) for the argument details.
Expand All @@ -175,29 +170,19 @@ https://en.wikipedia.org/wiki/Diversity_index#Shannon_index
function shannon_diversity(solution; threshold = eps(), kwargs...)
measure_on = extract_last_timesteps(solution; kwargs...)

if sum(measure_on) == 0
NaN
end
shan = [
shannon_diversity(vec(measure_on[:, i]); threshold = threshold) for
i in 1:size(measure_on, 2)
]
shan = shannon_diversity.(eachcol(measure_on), threshold = threshold)
mean(shan)

end

function shannon_diversity(n::AbstractVector; threshold = eps())
x = copy(n)
x = filter((k) -> k > threshold, x)
try
if length(x) >= 1
p = x ./ sum(x)
corr = log.(length(x))
p_ln_p = p .* log.(p)
-(sum(p_ln_p))
else
NaN
end
catch
x = filter(k -> k > threshold, n)

if length(x) >= 1
p = x ./ sum(x)
p_ln_p = p .* log.(p)
-(sum(p_ln_p))
else
NaN
end
end
Expand All @@ -206,7 +191,8 @@ end
"""
simpson(solution; threshold, kwargs...)
Computes the Simpson diversity index, i.e. the second Hill number.
Computes the average Simpson diversity index, i.e. the second Hill number,
over the `last` timesteps.
kwargs... arguments are dispatched to [`extract_last_timesteps`](@ref). See
[`extract_last_timesteps`](@ref) for the argument details.
Expand All @@ -219,27 +205,19 @@ https://en.wikipedia.org/wiki/Diversity_index#Simpson_index
"""
function simpson(solution; threshold = eps(), kwargs...)
measure_on = extract_last_timesteps(solution; kwargs...)
if sum(measure_on) == 0
NaN
end
simp = [
simpson(vec(measure_on[:, i]); threshold = threshold) for i in 1:size(measure_on, 2)
]

simp = simpson.(eachcol(measure_on), threshold = threshold)
mean(simp)
end

function simpson(n::AbstractVector; threshold = eps())
x = copy(n)
x = filter((k) -> k > threshold, x)
try
if length(x) >= 1
p = x ./ sum(x)
p2 = 2 .^ p
1 / sum(p2)
else
NaN
end
catch
x = filter((k) -> k > threshold, n)

if length(x) >= 1
p = x ./ sum(x)
p2 = 2 .^ p
1 / sum(p2)
else
NaN
end
end
Expand All @@ -248,7 +226,7 @@ end
"""
evenness(solution; threshold, kwargs...)
Computes the Pielou evenness over the `last` timesteps.
Computes the average Pielou evenness, over the `last` timesteps.
kwargs... arguments are dispatched to [`extract_last_timesteps`](@ref). See
[`extract_last_timesteps`](@ref) for the argument details.
Expand All @@ -261,25 +239,16 @@ https://en.wikipedia.org/wiki/Species_evenness
"""
function evenness(solution; threshold = eps(), kwargs...)
measure_on = extract_last_timesteps(solution; kwargs...)
if sum(measure_on) == 0
NaN
end
piel = [
evenness(vec(measure_on[:, i]); threshold = threshold) for
i in 1:size(measure_on, 2)
]

piel = evenness.(eachcol(measure_on), threshold = threshold)
mean(piel)
end

function evenness(n::AbstractVector; threshold = eps())
x = filter((k) -> k > threshold, n)
try
if length(x) > 0
shannon_diversity(x) / log(length(x))
else
NaN
end
catch
if length(x) > 0
shannon_diversity(x) / log(length(x))
else
NaN
end
end
Expand Down Expand Up @@ -406,7 +375,7 @@ julia> foodweb = FoodWeb([0 0 0; 0 1 0; 1 1 0]; quiet = true);
params = ModelParameters(foodweb);
B0 = [0.5, 0.5, 0.5];
sol = simulate(params, B0; verbose = false);
sum(trophic_structure(sol; last = 1).living_net) - sum(foodweb.A)
sum(trophic_structure(sol; last = 1).A) - sum(foodweb.A)
-2
```
"""
Expand All @@ -416,22 +385,32 @@ function trophic_structure(solution; threshold = eps(), idxs = nothing, kwargs..
`trophic_structure()` and it may never be. \
Please set `idxs = nothing`."))


living_sp = living_species(solution; threshold = threshold, kwargs...)

bm_sp = biomass(solution; kwargs...).sp[living_sp.idxs]

living_net = get_parameters(solution).network.A[living_sp.idxs, living_sp.idxs]
tlvl = trophic_levels(living_net)

(
max = maximum(tlvl),
avg = mean(tlvl),
w_avg = sum(tlvl .* (bm_sp ./ sum(bm_sp))),
living_species = living_sp.species,
tlvl = tlvl,
living_net = living_net,
)
net = get_parameters(solution).network.A[living_sp.idxs, living_sp.idxs]
tlvl = trophic_levels(net)

if isempty(tlvl)
(
max = NaN,
avg = NaN,
w_avg = NaN,
species = living_sp.species,
tlvl = tlvl,
A = net,
)
else
(
max = maximum(tlvl),
avg = mean(tlvl),
w_avg = sum(tlvl .* (bm_sp ./ sum(bm_sp))),
species = living_sp.species,
tlvl = tlvl,
A = net,
)
end
end

"""
Expand Down
92 changes: 70 additions & 22 deletions test/measures/test-functioning.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,15 +28,28 @@ end
@testset "Trophic structure" begin
foodweb = FoodWeb([0 1 1; 0 0 0; 0 0 0])
params = ModelParameters(foodweb)
sol = simulates(params, [0, 0, 0]; verbose = true)
sim_zero = simulates(params, [0, 0, 0]; verbose = true)
sim_three = simulates(params, [0.5, 0.5, 0.5]; verbose = true)

@test_throws ArgumentError("`idxs` is not currently supported for \
`trophic_structure()` and it may never be. \
Please set `idxs = nothing`.") trophic_structure(
sol;
sim_zero;
last = 10,
idxs = 2,
)

troph_zero = trophic_structure(sim_zero)
troph_three = trophic_structure(sim_three)
troph_three_nan = trophic_structure(sim_three; threshold = 1000)

@test [isnan(i) for i in troph_zero[(:max, :avg, :w_avg)]] |> all
@test [isempty(i) for i in troph_zero[(:species, :A, :tlvl)]] |> all
@test [isnan(i) for i in troph_three_nan[(:max, :avg, :w_avg)]] |> all
@test [isempty(i) for i in troph_three_nan[(:species, :A, :tlvl)]] |> all

@test troph_three[(:species, :A, :tlvl)] ==
(species = foodweb.species, A = foodweb.A, tlvl = [2.0, 1.0, 1.0])
end

@testset "Producer growth rate" begin
Expand Down Expand Up @@ -87,41 +100,49 @@ end
sim_zero = simulates(params, [0, 0]; verbose = false)

# Total biomass should converge to K
@test biomass(sim_one_sp; last = 1).total ≈ get_parameters(sim_one_sp).environment.K[2] rtol =
0.001
tmp_K = get_parameters(sim_one_sp).environment.K[2]
@test biomass(sim_one_sp; last = 1).total ≈ tmp_K rtol = 0.001
@test biomass(sim_one_sp; last = 1).total ≈ tmp_K rtol = 0.001
bm_two_sp = biomass(sim_two_sp; last = 1)
@test bm_two_sp.total == sum(bm_two_sp.sp)
@test bm_two_sp.total ≈ 2 atol = 0.001

# Species sub selection works
@test biomass(sim_one_sp; last = 1, idxs = [1]).sp[1] ≈
biomass(sim_one_sp; last = 1, idxs = [1]).total

@test biomass(sim_zero; last = 2).total ≈ 0.0
# Species richness is the binary equivalent of total_biomass
@test EcologicalNetworksDynamics.richness(sim_zero; last = 2) ≈ 0.0
@test richness(sim_zero; last = 2) ≈ 0.0
@test species_persistence(sim_zero; last = 2) ≈ 0.0 # Weird but it is a feature

@test EcologicalNetworksDynamics.richness(sim_two_sp[:, end]) ==
EcologicalNetworksDynamics.richness(sim_two_sp; last = 1) ==
2
@test richness(sim_two_sp[:, end]) == richness(sim_two_sp; last = 1) == 2

@test EcologicalNetworksDynamics.richness(sim_one_sp[:, end]) ==
EcologicalNetworksDynamics.richness(sim_one_sp; last = 1) ==
1
@test richness(sim_one_sp[:, end]) == richness(sim_one_sp; last = 1) == 1

@test EcologicalNetworksDynamics.richness(sim_zero[:, end]) ==
EcologicalNetworksDynamics.richness(sim_zero; last = 1) ==
0
@test richness(sim_zero[:, end]) == richness(sim_zero; last = 1) == 0

# Other hill diversity numbers
## Shannon
@test shannon_diversity(sim_two_sp[:, end]) ==
shannon_diversity(sim_two_sp; last = 1) ==
shannon_diversity([1, 1]) ==
log(2)

# 0 entropy if 1 species only
@test shannon_diversity(sim_one_sp[:, end]) ==
shannon_diversity(sim_one_sp; last = 1) ≈
0.0 # 0 entropy if 1 species only

shannon_diversity([1]) ==
shannon_diversity(sim_two_sp; idxs = 1) ==
shannon_diversity(sim_two_sp; idxs = "s1") ==
shannon_diversity(sim_two_sp; idxs = 2) ==
shannon_diversity(sim_two_sp; idxs = "s2") ==
shannon_diversity(sim_one_sp; last = 1) ==
0.0

# Not defined for 0 species
@test isnan(shannon_diversity(sim_zero[:, end])) ==
isnan(shannon_diversity([0, 0])) ==
isnan(shannon_diversity([1, 1]; threshold = 1)) ==
isnan(shannon_diversity(sim_zero; last = 1)) # Not defined for 0 species

@test shannon_diversity(sim_one_sp[:, end]) ==
Expand All @@ -134,16 +155,43 @@ end
## Simpson
@test simpson(sim_two_sp[:, end]) ==
simpson(sim_two_sp; last = 1) ==
simpson([1, 1]) ==
1 / sum(2 .^ [1 / 2, 1 / 2])
# 0.5 if 1 species only
@test simpson(sim_one_sp[:, end]) ==
simpson(sim_one_sp; last = 1) ==
simpson(sim_two_sp; idxs = 1) ==
simpson(sim_two_sp; idxs = "s1") ==
simpson(sim_two_sp; idxs = 2) ==
simpson(sim_two_sp; idxs = "s2") ==
simpson([1]) ==
1 / sum(2 .^ 1) ==
0.5# 0.5 if 1 species only
@test isnan(simpson(sim_zero[:, end])) == isnan(simpson(sim_zero; last = 1))
0.5
@test isnan(simpson(sim_zero[:, end])) ==
isnan(simpson(sim_zero; last = 1)) ==
isnan(simpson([0, 0])) ==
isnan(simpson([1, 1]; threshold = 1))

# Community evenness
@test evenness(sim_two_sp[:, end]) == evenness(sim_two_sp; last = 1) ≈ 1.0 # Maximum equity of species biomass
@test isnan(evenness(sim_one_sp[:, end])) == isnan(evenness(sim_one_sp; last = 1))# Should be NaN if 1 species
@test isnan(evenness(sim_zero[:, end])) == isnan(evenness(sim_zero; last = 1))# Should be NaN if 0 species
# Maximum equity of species biomass
@test evenness(sim_two_sp[:, end]) ==
evenness(sim_two_sp; last = 1) ==
evenness([1, 1]) ==
1.0
# Should be NaN if 1 species
@test isnan(evenness(sim_one_sp[:, end])) ==
isnan(evenness(sim_one_sp; last = 1)) ==
isnan(evenness(sim_two_sp; idxs = 1)) ==
isnan(evenness(sim_two_sp; idxs = "s1")) ==
isnan(evenness(sim_two_sp; idxs = 2)) ==
isnan(evenness(sim_two_sp; idxs = "s2")) ==
isnan(evenness([1, 0])) ==
isnan(evenness([2, 1]; threshold = 1))
# Should be NaN if 0 species
@test isnan(evenness(sim_zero[:, end])) ==
isnan(evenness(sim_zero; last = 1)) ==
isnan(evenness([0, 0])) ==
isnan(evenness([1, 1]; threshold = 1))


end
3 changes: 1 addition & 2 deletions test/measures/test-utils.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,14 +51,13 @@
)

@test_throws ArgumentError("Cannot extract timesteps with \
Matrix{Any}(undef, 0, 0). \
Any[]. \
`last` should be an Integer or a String \
representing a percentage.") extract_last_timesteps(
sim,
last = [],
)


msg_warn = "Cannot extract 100 timesteps from a trajectory solution with only 12 \
timesteps. We set 'last' to 1, i.e. the last 10% of timesteps."
@test_warn msg_warn extract_last_timesteps(sim, last = 100, autofix_last = true)
Expand Down

0 comments on commit 07e68a9

Please sign in to comment.