From e079ff51c0d3d681b4ce21c10a30f8ff2a4f2ea3 Mon Sep 17 00:00:00 2001
From: Alain Danet <alain.danet@mnhn.fr>
Date: Sat, 4 Feb 2023 12:32:50 +0000
Subject: [PATCH] Improve tests and output utils according to Iago review
 :tornado::

- Reference: https://github.com/BecksLab/BEFWM2/pull/87#pullrequestreview-1272831620
---
 src/EcologicalNetworksDynamics.jl |   3 +
 src/measures/functioning.jl       | 237 +++++++++++++++---------------
 src/measures/stability.jl         | 129 ++++++++++------
 src/measures/utils.jl             | 126 +++++++++++-----
 test/measures/test-functioning.jl |  43 +++++-
 test/measures/test-stability.jl   |  65 ++++++--
 test/measures/test-utils.jl       |  90 ++++++++++--
 test/runtests.jl                  |   3 +-
 8 files changed, 462 insertions(+), 234 deletions(-)

diff --git a/src/EcologicalNetworksDynamics.jl b/src/EcologicalNetworksDynamics.jl
index 558fc0828..24b32e88f 100644
--- a/src/EcologicalNetworksDynamics.jl
+++ b/src/EcologicalNetworksDynamics.jl
@@ -63,6 +63,8 @@ export ClassicResponse
 export connectance
 export cpad
 export coefficient_of_variation
+export community_cv
+export cv_sp
 export DefaultGrowthParams
 export DefaultMaxConsumptionParams
 export DefaultMetabolismParams
@@ -126,6 +128,7 @@ export richness
 export set_temperature!
 export shannon_diversity
 export simpson
+export synchrony
 export simulate
 export species_persistence
 export trophic_levels
diff --git a/src/measures/functioning.jl b/src/measures/functioning.jl
index f7ba5d4ec..86d9f4952 100644
--- a/src/measures/functioning.jl
+++ b/src/measures/functioning.jl
@@ -8,7 +8,7 @@ Adapted from BioenergeticFoodWeb.jl
 
 Returns the average number of species with a biomass larger than `threshold`
 over the `last` timesteps. `kwargs...` are optional arguments passed to
-`extract_last_timesteps`.
+[`extract_last_timesteps`](@ref).
 
 # Arguments:
 
@@ -19,14 +19,10 @@ over the `last` timesteps. `kwargs...` are optional arguments passed to
 
 ```jldoctest
 julia> foodweb = FoodWeb([0 0; 1 0]);
-
-julia> params = ModelParameters(foodweb);
-
-julia> B0 = [0.5, 0.5];
-
-julia> sol = simulate(params, B0);
-
-julia> richness(sol; last = 10)
+       params = ModelParameters(foodweb);
+       B0 = [0.5, 0.5];
+       sol = simulate(params, B0);
+       richness(sol; last = 10)
 2.0
 
 julia> sha = shannon_diversity(sol);
@@ -46,14 +42,14 @@ function richness(solution; threshold = eps(), kwargs...)
     measure_on = extract_last_timesteps(solution; kwargs...)
 
     rich = [
-        BEFWM2.richness(vec(measure_on[:, i]); threshold = threshold) for
+        richness(vec(measure_on[:, i]); threshold = threshold) for
         i in 1:size(measure_on, 2)
     ]
     mean(rich)
 end
 
 """
-    richness(n::Vector; threshold = eps())
+    richness(n::AbstractVector; threshold = eps())
 
 When applied to a vector of biomass, returns the number of biomass above `threshold`
 
@@ -67,18 +63,16 @@ julia> richness([1, 1])
 2
 ```
 """
-richness(n::Vector; threshold = eps()) = sum(n .> threshold)
+richness(n::AbstractVector; threshold = eps()) = sum(n .> threshold)
 
 """
     species_persistence(solution; threshold, kwargs...)
 
-Returns the average proportion of living species over the `last` timesteps.
+Returns the average proportion of species having a biomass superior or equal to threshold
+over the `last` timesteps.
 
-See [`richness`](@ref) for the argument details.
-
-When applied to a vector of biomass, e.g.
-`species_persistence(n::Vector; threshold = eps())`, it returns the proportion of
-species which biomass is above `threshold`.
+kwargs... arguments are dispatched to [`extract_last_timesteps`](@ref). See
+[`extract_last_timesteps`](@ref) for the argument details.
 
 When applied to a vector of biomass, e.g.
 `species_persistence(n::Vector; threshold = eps())`, it returns the proportion of
@@ -87,23 +81,18 @@ species which biomass is above `threshold`.
 # Examples
 
 ```jldoctest
-julia> foodweb = FoodWeb([0 0; 0 0]; quiet = true); # A foodweb with two producers.
-
-julia> params = ModelParameters(foodweb);
-
-julia> sim_two = simulate(params, [0.5, 0.5]);
-
-julia> species_persistence(sim_two; last = 1) # Both producers survived.
+julia> foodweb = FoodWeb([0 0; 0 0]; quiet = true);
+       params = ModelParameters(foodweb);
+       sim_two = simulate(params, [0.5, 0.5]);
+       species_persistence(sim_two; last = 1)
 1.0
 
 julia> sim_one = simulate(params, [0, 0.5]);
-
-julia> species_persistence(sim_one; last = 1) # Only one producer survived.
+       species_persistence(sim_one; last = 1)
 0.5
 
 julia> sim_zero = simulate(params, [0, 0]);
-
-julia> species_persistence(sim_zero; last = 1) # Degenerated situation.
+       species_persistence(sim_zero; last = 1)
 0.0
 
 julia> species_persistence([0, 1])
@@ -118,7 +107,7 @@ function species_persistence(solution; kwargs...)
     m = size(solution, 1) # Number of species is the number of rows in the biomass matrix
     r / m
 end
-species_persistence(n::Vector; threshold = eps()) =
+species_persistence(n::AbstractVector; threshold = eps()) =
     richness(n; threshold = threshold) / length(n)
 
 """
@@ -129,10 +118,12 @@ Returns a named tuple of total and species biomass, averaged over the last `last
 # Arguments
 
   - `solution`, `last`: See [`richness`](@ref) doc.
-  - `idxs`: See [`extract_last_timesteps`](@ref) doc.
 
-Can also handle a species x time biomass matrix, e.g. `biomass(mat::Matrix;)` or a
-vector, e.g. `biomass(vec::Vector;)`.
+kwargs... arguments are dispatched to [`extract_last_timesteps`](@ref). See
+[`extract_last_timesteps`](@ref) for the argument details.
+
+Can also handle a species x time biomass matrix, e.g. `biomass(mat::AbstractMatrix;)` or a
+vector, e.g. `biomass(vec::AbstractVector;)`.
 
 Can also handle a species x time biomass matrix, e.g. `biomass(mat::Matrix;)` or a
 vector, e.g. `biomass(vec::Vector;)`.
@@ -141,14 +132,10 @@ vector, e.g. `biomass(vec::Vector;)`.
 
 ```jldoctest
 julia> foodweb = FoodWeb([0 0; 1 0]);
-
-julia> params = ModelParameters(foodweb);
-
-julia> B0 = [0.5, 0.5];
-
-julia> sol = simulate(params, B0);
-
-julia> bm = biomass(sol; last = 2);
+       params = ModelParameters(foodweb);
+       B0 = [0.5, 0.5];
+       sol = simulate(params, B0);
+       bm = biomass(sol; last = 2);
        biomass(sol; last = 2).sp ≈ [0.1890006203352691, 0.21964742227673806]
 true
 
@@ -167,17 +154,19 @@ function biomass(solution; kwargs...)
     measure_on = extract_last_timesteps(solution; kwargs...)
     biomass(measure_on)
 end
-biomass(mat::Matrix;) = (total = mean(vec(sum(mat; dims = 1))), sp = mean.(eachrow(mat)))
-biomass(vec::Vector;) = (total = mean(vec), sp = mean(vec))
+biomass(mat::AbstractMatrix;) =
+    (total = mean(vec(sum(mat; dims = 1))), sp = mean.(eachrow(mat)))
+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.
 
-See [`richness`](@ref) doc.
+kwargs... arguments are dispatched to [`extract_last_timesteps`](@ref). See
+[`extract_last_timesteps`](@ref) for the argument details.
 
-Can also handle a vector, e.g. shannon_diversity(n::Vector; threshold = eps())
+Can also handle a vector, e.g. shannon_diversity(n::AbstractVector; threshold = eps())
 
 # Reference
 
@@ -196,7 +185,7 @@ function shannon_diversity(solution; threshold = eps(), kwargs...)
     mean(shan)
 end
 
-function shannon_diversity(n::Vector; threshold = eps())
+function shannon_diversity(n::AbstractVector; threshold = eps())
     x = copy(n)
     x = filter((k) -> k > threshold, x)
     try
@@ -219,16 +208,16 @@ end
 
 Computes the Simpson diversity index, i.e. the second Hill number.
 
-See [`richness`](@ref) doc.
+kwargs... arguments are dispatched to [`extract_last_timesteps`](@ref). See
+[`extract_last_timesteps`](@ref) for the argument details.
 
-Can also handle a vector, e.g. simpson(n::Vector; threshold = eps())
+Can also handle a vector, e.g. simpson(n::AbstractVector; threshold = eps())
 
 # Reference
 
 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
@@ -239,7 +228,7 @@ function simpson(solution; threshold = eps(), kwargs...)
     mean(simp)
 end
 
-function simpson(n::Vector; threshold = eps())
+function simpson(n::AbstractVector; threshold = eps())
     x = copy(n)
     x = filter((k) -> k > threshold, x)
     try
@@ -261,9 +250,10 @@ end
 
 Computes the Pielou evenness over the `last` timesteps.
 
-See [`richness`](@ref) doc.
+kwargs... arguments are dispatched to [`extract_last_timesteps`](@ref). See
+[`extract_last_timesteps`](@ref) for the argument details.
 
-Can also handle a vector, e.g. evenness(n::Vector; threshold = eps())
+Can also handle a vector, e.g. `evenness(n::AbstractVector; threshold = eps())`
 
 # Reference
 
@@ -281,7 +271,7 @@ function evenness(solution; threshold = eps(), kwargs...)
     mean(piel)
 end
 
-function evenness(n::Vector; threshold = eps())
+function evenness(n::AbstractVector; threshold = eps())
     x = filter((k) -> k > threshold, n)
     try
         if length(x) > 0
@@ -300,22 +290,18 @@ end
 Returns growth rates of producers over the `last` timesteps, as well as the
 average and the standard deviation.
 
-See [`richness`](@ref) doc for argument description.
+kwargs... arguments are dispatched to [`extract_last_timesteps`](@ref). See
+[`extract_last_timesteps`](@ref) for the argument details.
 
 # Examples
 
 ```jldoctest
 julia> foodweb = FoodWeb([0 1 1; 0 0 0; 0 0 0]);
-
-julia> params = ModelParameters(foodweb);
-
-julia> B0 = [0.5, 0.5, 0.5];
-
-julia> sol = simulate(params, B0);
-
-julia> g = producer_growth(sol; last = 10);
-
-julia> g.species, round.(g.mean, digits = 2)
+       params = ModelParameters(foodweb);
+       B0 = [0.5, 0.5, 0.5];
+       sol = simulate(params, B0);
+       g = producer_growth(sol; last = 10);
+       g.species, round.(g.mean, digits = 2)
 (["s2", "s3"], [0.15, 0.15])
 ```
 """
@@ -328,10 +314,10 @@ function producer_growth(solution; kwargs...)
 
     Kp = parameters.environment.K[producer_idxs]
     rp = parameters.biorates.r[producer_idxs]
-    αp = parameters.producer_competition.α[producer_idxs, producer_idxs] # matrix of producer competition
+    αp = parameters.producer_competition.α[producer_idxs, producer_idxs]
 
     #Extract the producer_species biomass over the last timesteps
-    measure_on = extract_last_timesteps(solution; kwargs...)[producer_idxs, :]
+    measure_on = extract_last_timesteps(solution; idxs = producer_idxs, kwargs...)
 
     # Compute growth and make it a time x species matrix
     #
@@ -362,12 +348,9 @@ Extract list of extinct species from the solution returned by `simulate()`.
 
 ```jldoctest
 julia> foodweb = FoodWeb([0 0; 1 0]);
-
-julia> params = ModelParameters(foodweb);
-
-julia> sol = simulate(params, [0.5, 0.5]);
-
-julia> isempty(get_extinct_species(sol)) # no species extinct
+       params = ModelParameters(foodweb);
+       sol = simulate(params, [0.5, 0.5]);
+       isempty(get_extinct_species(sol)) # no species extinct
 true
 ```
 
@@ -382,12 +365,9 @@ Extract the [`ModelParameters`](@ref) input from the solution returned by `simul
 
 ```jldoctest
 julia> foodweb = FoodWeb([0 0; 1 0]);
-
-julia> params = ModelParameters(foodweb);
-
-julia> sol = simulate(params, [0.5, 0.5]);
-
-julia> isa(get_parameters(sol), ModelParameters)
+       params = ModelParameters(foodweb);
+       sol = simulate(params, [0.5, 0.5]);
+       isa(get_parameters(sol), ModelParameters)
 true
 ```
 
@@ -396,85 +376,82 @@ See also [`simulate`](@ref), [`get_extinct_species`](@ref).
 get_parameters(sol) = sol.prob.p.params
 
 """
-    trophic_structure(solution; threshold, last, idxs = nothing)
+    trophic_structure(solution; threshold, kwargs...)
 
 Returns the maximum trophic level, its average and its biomass weighted average,
 computed over the living species (see [`living_species`](@ref)).  It also returns the
 adjacency matrix containing only the living species and the vector of the living species.
 
-See [`richness`](@ref) for the argument details.
+kwargs... arguments are dispatched to [`extract_last_timesteps`](@ref). See
+[`extract_last_timesteps`](@ref) for the argument details.
+
+# Examples
 
 ```jldoctest
 julia> foodweb = FoodWeb([0 1 1; 0 0 0; 0 0 0]);
-
-julia> params = ModelParameters(foodweb);
-
-julia> B0 = [0.5, 0.5, 0.5];
-
-julia> sol = simulate(params, B0; verbose = true);
-
-julia> three_sp = trophic_structure(sol; last = 10);
-
-julia> three_sp[(:max, :avg)]
+       params = ModelParameters(foodweb);
+       B0 = [0.5, 0.5, 0.5];
+       sol = simulate(params, B0; verbose = true);
+       three_sp = trophic_structure(sol; last = 10);
+       three_sp[(:max, :avg)]
 (max = 2.0, avg = 1.3333333333333333)
 
 julia> B0 = [0, 0.5, 0.5];
-
-julia> sol = simulate(params, B0; verbose = true);
-
-julia> no_consumer = trophic_structure(sol; last = 10);
-
-julia> no_consumer[(:max, :avg)]
+       sol = simulate(params, B0; verbose = true);
+       no_consumer = trophic_structure(sol; last = 10);
+       no_consumer[(:max, :avg)]
 (max = 1.0, avg = 1.0)
 
 julia> foodweb = FoodWeb([0 0 0; 0 1 0; 1 1 0]; quiet = true);
-
-julia> params = ModelParameters(foodweb);
-
-julia> B0 = [0.5, 0.5, 0.5];
-
-julia> sol = simulate(params, B0; verbose = false);
-
-julia> sum(trophic_structure(sol; last = 1).living_net) - sum(foodweb.A) #2 links lost during simulation
+       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)
 -2
 ```
 """
-function trophic_structure(solution; threshold = eps(), kwargs...)
+function trophic_structure(solution; threshold = eps(), idxs = nothing, kwargs...)
+
+    isnothing(idxs) || throw(ArgumentError("`idxs` is not currently supported for \
+                                           `trophic_structure()` and it may never be. \
+                                           Please set `idxs = nothing`."))
+
 
-    bm_sp = biomass(solution; kwargs...).sp
     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[living_sp.idxs] ./ sum(bm_sp[living_sp.idxs]))),
+        w_avg = sum(tlvl .* (bm_sp ./ sum(bm_sp))),
         living_species = living_sp.species,
+        tlvl = tlvl,
         living_net = living_net,
     )
 end
 
 """
-    living_species(solution; threshold = eps(), kwargs...)
+    living_species(solution; threshold, kwargs...)
 
 Returns the vectors of alive species and their indices in the original network.
 Living species are the ones having, in average, a biomass above `threshold` over
-the `last` timesteps.
+the `last` timesteps. `kwargs...` are optional arguments passed to
+[`extract_last_timesteps`](@ref).
 
-See [`richness`](@ref) for the argument details.
+kwargs... arguments are dispatched to [`extract_last_timesteps`](@ref). See
+[`extract_last_timesteps`](@ref) for the argument details.
 
 # Examples
 
 ```jldoctest
 julia> foodweb = FoodWeb([0 1 1; 0 0 0; 0 0 0]);
-
-julia> params = ModelParameters(foodweb);
-
-julia> B0 = [0.5, 0.5, 0.5];
-
-julia> sol = simulate(params, B0; verbose = true);
+       params = ModelParameters(foodweb);
+       B0 = [0.5, 0.5, 0.5];
+       sol = simulate(params, B0; verbose = true);
 
 julia> living_species(sol)
 (species = ["s1", "s2", "s3"], idxs = [1, 2, 3])
@@ -487,20 +464,36 @@ julia> living_species(sol; last = 1)
 (species = ["s2", "s3"], idxs = [2, 3])
 ```
 """
-function living_species(solution; threshold = eps(), kwargs...)
+function living_species(solution; threshold = eps(), idxs = nothing, kwargs...)
+
+    measure_on = extract_last_timesteps(solution; idxs = idxs, kwargs...)
 
-    bm_sp = biomass(solution; kwargs...).sp
-    living_sp = findall(x -> x >= threshold, bm_sp)
+    alive_sp = living_species(measure_on; threshold = threshold)
 
-    (species = get_parameters(solution).network.species[living_sp], idxs = living_sp)
+    sp = get_parameters(solution).network.species
+
+    if !isnothing(idxs) # if idxs has been used in extract_last_timesteps
+        idxs = process_idxs(solution, idxs)
+        alive_sp = idxs[alive_sp]
+        sp_alive_name = sp[alive_sp]
+    else
+        sp_alive_name = sp[alive_sp]
+    end
+
+    (species = sp_alive_name, idxs = alive_sp)
 end
 
+living_species(mat::Matrix; threshold = eps()) =
+    findall(x -> x >= threshold, biomass(mat).sp)
+
 """
     min_max(solution; kwargs...)
 
-Returns the vectors of minimum and maximum biomass of each species over the `last` timesteps.
+Returns the vectors of minimum and maximum biomass of each species over the `last`
+timesteps.
 
-See [`richness`](@ref) for the argument details.
+kwargs... arguments are dispatched to [`extract_last_timesteps`](@ref). See
+[`extract_last_timesteps`](@ref) for the argument details.
 
 # Examples
 
diff --git a/src/measures/stability.jl b/src/measures/stability.jl
index fa52bf9c5..11f1feb89 100644
--- a/src/measures/stability.jl
+++ b/src/measures/stability.jl
@@ -3,38 +3,51 @@ Various measures of stability
 =#
 
 """
-    cv_sp(mat)
+    cv_sp(solution; threshold, kwargs...)
 
 Computes the temporal CV of species biomass and its average.
 
-Computes the temporal CV of species biomass and its average.
+See [`coefficient_of_variation`](@ref) for the details.
+"""
+function cv_sp(solution; threshold = eps(), kwargs...)
 
-  - mat: A time x species biomass matrix. Typically the transposition of the
-    output of `DifferentialEquations.solve()` or `BEFWM2.simulate()`)
+    measure_on = extract_last_timesteps(solution; kwargs...)
+    # Fetch species that are alive, whose avg biomass is > threshold
+    living_sp = living_species(measure_on; threshold = threshold)
 
-See [`cv`](@ref) for the details.
-"""
-function cv_sp(mat)
+    # Transpose to get the time x species matrix
+    mat = transpose(measure_on[living_sp, :])
+
+    cvsp = cv_sp(mat)
+
+    out = (avg_cv_sp = cvsp.avg, cv_sp = cvsp.sp)
 
-    avg_sp = mean.(eachcol(mat))
-    std_sp = std.(eachcol(mat), corrected = false)
+    out
+end
+function cv_sp(mat::AbstractMatrix)
 
-    rel_sp = avg_sp ./ sum(avg_sp)
-    rel_sd_sp = std_sp ./ avg_sp
+    if any(size(mat) .== 0)
+        cv = (avg = NaN, sp = NaN)
+    else
+        avg_sp = mean.(eachcol(mat))
+        std_sp = std.(eachcol(mat), corrected = false)
 
-    avg_cv_sp = sum(rel_sp .* rel_sd_sp)
+        rel_sp = avg_sp ./ sum(avg_sp)
+        rel_sd_sp = std_sp ./ avg_sp
 
-    (avg = avg_cv_sp, sp = std_sp ./ avg_sp)
+        avg_cv_sp = sum(rel_sp .* rel_sd_sp)
 
+        cv = (avg = avg_cv_sp, sp = std_sp ./ avg_sp)
+    end
 end
 
 """
-    synchrony(mat)
+    synchrony(solution; threshold, kwargs...)
 
 Compute the synchrony of species biomass fluctuations following
 Loreau & de Mazancourt (2008).
 
-See [`avg_cv_sp`](@ref) for the argument details.
+See [`coefficient_of_variation`](@ref) for the argument details.
 
 # Reference
 
@@ -42,37 +55,69 @@ Loreau, M., & de Mazancourt, C. (2008). Species Synchrony and Its Drivers :
 Neutral and Nonneutral Community Dynamics in Fluctuating Environments. The
 American Naturalist, 172(2), E48‑E66. https://doi.org/10.1086/589746
 """
-function synchrony(mat)
+function synchrony(solution; threshold = eps(), kwargs...)
+
+    measure_on = extract_last_timesteps(solution; kwargs...)
+    # Fetch species that are alive, whose avg biomass is > threshold
+    living_sp = living_species(measure_on; threshold = threshold)
+
+    # Transpose to get the time x species matrix
+    mat = transpose(measure_on[living_sp, :])
+
+    synchrony(mat)
+
+end
+function synchrony(mat::AbstractMatrix)
 
-    cov_mat = cov(mat; corrected = false)
+    if any(size(mat) .== 0)
+        phi = NaN
+    else
+        cov_mat = cov(mat; corrected = false)
 
-    com_var = sum(cov_mat)
-    std_sp = sum(std.(eachcol(mat), corrected = false))
+        com_var = sum(cov_mat)
+        std_sp = sum(std.(eachcol(mat), corrected = false))
 
-    phi = com_var / std_sp^2
+        phi = com_var / std_sp^2
+    end
 
     phi
+
 end
 
 """
-    temporal_cv(mat)
+    community_cv(solution; threshold, kwargs...)
 
-Compute the temporal CV of community biomass.
+Compute the temporal Coefficient of Variation of community biomass.
 
-See [`avg_cv_sp`](@ref) for the argument details.
+See [`coefficient_of_variation`](@ref) for the argument details.
 """
-function temporal_cv(mat)
+function community_cv(solution; threshold = eps(), kwargs...)
+
+    measure_on = extract_last_timesteps(solution; kwargs...)
+    # Fetch species that are alive, whose avg biomass is > threshold
+    living_sp = living_species(measure_on; threshold = threshold, kwargs...)
+
+    # Transpose to get the time x species matrix
+    mat = transpose(measure_on[living_sp, :])
 
-    total_com_bm = sum.(eachrow(mat))
-    cv_com = std(total_com_bm; corrected = false) / mean(total_com_bm)
+    community_cv(mat)
 
+end
+function community_cv(mat::AbstractMatrix)
+
+    if any(size(mat) .== 0)
+        cv_com = NaN
+    else
+        total_com_bm = sum.(eachrow(mat))
+        cv_com = std(total_com_bm; corrected = false) / mean(total_com_bm)
+    end
     cv_com
 end
 
 """
     coefficient_of_variation(solution; threshold, last)
 
-Computes Coefficient of Variation (CV) of community biomass and its partition in
+Computes the Coefficient of Variation (CV) of community biomass and its partition in
 average species CV and synchrony, following Thibault & Connolly (2013).
 
 The function excludes dead species, i.e. species that have an average biomass
@@ -91,25 +136,17 @@ relationships : Towards a unified model of portfolio effects. Ecology Letters,
 
 ```jldoctest
 julia> foodweb = FoodWeb([0 1 1; 0 0 0; 0 0 0]); # Two producers and one consumer
-
-julia> params = ModelParameters(foodweb);
-
-julia> B0 = [0.5, 0.5, 0.5];
-
-julia> sol = simulate(params, B0; verbose = true);
-
-julia> s = coefficient_of_variation(sol; last = 10);
-
-julia> map(x -> round(x; digits = 2), s[(:cv_com, :avg_cv_sp, :synchrony)])
+       params = ModelParameters(foodweb);
+       B0 = [0.5, 0.5, 0.5];
+       sol = simulate(params, B0; verbose = true);
+       s = coefficient_of_variation(sol; last = 10);
+       map(x -> round(x; digits = 2), s[(:cv_com, :avg_cv_sp, :synchrony)])
 (cv_com = 0.02, avg_cv_sp = 0.06, synchrony = 0.1)
 
 julia> B0 = [0, 0.5, 0.5]; # Two producers
-
-julia> sol = simulate(params, B0; verbose = true);
-
-julia> s = coefficient_of_variation(sol; last = 10);
-
-julia> map(x -> round(x; digits = 2), s[(:cv_com, :avg_cv_sp, :synchrony)])
+       sol = simulate(params, B0; verbose = true);
+       s = coefficient_of_variation(sol; last = 10);
+       map(x -> round(x; digits = 2), s[(:cv_com, :avg_cv_sp, :synchrony)])
 (cv_com = 0.14, avg_cv_sp = 0.14, synchrony = 1.0)
 ```
 """
@@ -117,14 +154,14 @@ function coefficient_of_variation(solution; threshold = eps(), kwargs...)
 
     measure_on = extract_last_timesteps(solution; kwargs...)
     # Fetch species that are alive, whose avg biomass is > threshold
-    living_sp = living_species(solution; threshold = threshold, kwargs...)
+    alive_sp = living_species(measure_on; threshold = threshold)
 
     # Transpose to get the time x species matrix
-    mat = transpose(measure_on[living_sp.idxs, :])
+    mat = transpose(measure_on[alive_sp, :])
 
     cvsp = cv_sp(mat)
     sync = synchrony(mat)
-    cv_com = temporal_cv(mat)
+    cv_com = community_cv(mat)
 
     out = (cv_com = cv_com, avg_cv_sp = cvsp.avg, synchrony = sync, cv_sp = cvsp.sp)
 
diff --git a/src/measures/utils.jl b/src/measures/utils.jl
index c174eb86f..63e7e9d94 100644
--- a/src/measures/utils.jl
+++ b/src/measures/utils.jl
@@ -10,7 +10,7 @@ Returns the biomass matrix of species x time over the `last` timesteps.
 # Arguments
 
   - `last`: the number of last timesteps to consider. A percentage can also be also be
-    provided as a `Char` ending by `%`. Defaulted to 10%.
+    provided as a `String` ending by `%`. Defaulted to 10%.
   - `idxs`: vector of species indexes or names. Set to `nothing` by default.
   - `autofix_last`: fixes wrong `last` specification if `last` was specified as an integer.
   - `quiet`: ignores warning issued by `autofix_last`.
@@ -33,7 +33,7 @@ true
 julia> sim = extract_last_timesteps(m; last = 1, idxs = [1, 2]);
        sim ==
        extract_last_timesteps(m; last = 1, idxs = ["s1", "s2"]) ≈
-       [0.188980349; 0.219659439][:, [1]]
+       [0.188980349; 0.219659439;;]
 true
 
 julia> sim = extract_last_timesteps(m; last = 1, idxs = [2]);
@@ -45,47 +45,103 @@ julia> extract_last_timesteps(m; last = 1000, autofix_last = true, quiet = true)
 true
 ```
 """
-function extract_last_timesteps(
-    solution;
-    last = "10%",
-    idxs = nothing,
-    autofix_last = false,
-    quiet = false,
-)
+function extract_last_timesteps(solution; idxs = nothing, kwargs...)
+
+    last = process_last_timesteps(solution; kwargs...)
+
+    out = solution[:, end-(last-1):end]
+
+    if !isnothing(idxs)
+        idxs = process_idxs(solution, idxs;)
+        # Extract species indices:
+        out = out[idxs, :]
+    end
+    deepcopy(out)
+end
+
+"""
+process_idxs(solution, idxs;)
+
+Check and sanatize the species indices or names provided (`idxs`). Used in
+[`extract_last_timesteps`](@ref) and [`living_species`](@ref).
+"""
+function process_idxs(solution, idxs;)
+    # Handle species names:
+    sp = get_parameters(solution).network.species
+    if eltype(idxs) <: AbstractString || eltype(idxs) <: Char
+        if idxs isa String
+            idxs = [idxs]
+        end
+
+        idxs_in = indexin(idxs, sp)
+        # Handle missing species
+        absent_sp = [isnothing(i) for i in idxs_in]
+        !any(absent_sp) || throw(ArgumentError("Species $(idxs[absent_sp]) are not \
+                                               found in the network. Any mispelling?"))
+
+        idxs = idxs_in[.!absent_sp]
+
+    elseif eltype(idxs) <: Integer
+        if idxs isa Integer
+            idxs = [idxs]
+        end
+        check_bound_idxs = map(x -> x >= 1 && x <= length(sp), idxs)
+        absent_sp = idxs[findall(.!check_bound_idxs)]
+        all(check_bound_idxs) || throw(ArgumentError("Cannot extract idxs \
+                                                     $(absent_sp). Each species index \
+                                                     should  be superior or equal to 1 \
+                                                     and should be inferior or equal \
+                                                     to the number of species \
+                                                     ($(length(sp)))."))
+
+    else
+        throw(ArgumentError("`idxs` should be a vector of Integer (species indices) or \
+                            String (species names)"))
+    end
+
+    idxs
+end
+
+"""
+"""
+function process_last_timesteps(solution; last = "10%", autofix_last = false, quiet = false)
 
     n_timesteps = length(solution.t)
-    @assert any(eltype(last) .=== [Int, Char]) "`last` should be an integer or a character"
 
-    # If last is a percentage
-    if eltype(last) === Char
-        @assert occursin("%", last) "`last` Char should end by %"
+    if last isa AbstractString
+        endswith(last, "%") || throw(ArgumentError("`last` String should end by %."))
         perc = parse(Float64, last[1:(end-1)])
-        @assert 0 < perc <= 100 "`last` should be > 0 and <= 100"
+        is_valid_perc = perc > 0.0 && perc <= 100.0
+        is_valid_perc || throw(ArgumentError("Cannot extract $(perc)% of the solution's \
+                                             timesteps: 0% < `last` <= 100% must hold."))
         last = round(Int, n_timesteps * perc / 100)
+    elseif last isa Integer
+        last > 0 || throw(ArgumentError("Cannot extract $last timesteps. `last` should be \
+                                        an Integer strictly superior to 0."))
+    elseif last isa Float64
+        throw(ArgumentError("Cannot extract `last` from a decimal number. \
+                            Did you mean \"$last%\"?"))
+    else
+        throw(ArgumentError("Cannot extract timesteps with $last. \
+                            `last` should be an Integer or a String \
+                            representing a percentage."))
     end
 
     if last > n_timesteps
-        @assert autofix_last "`last` > length(solution.t),
-        please decrease `last` or consider to setting `autofix_last = true`."
-
-        last = round(Int, n_timesteps * 10 / 100)
-        quiet || @warn "'last' is lower or equal to the number of timesteps. \
-        'last' is then set to $(last), i.e. the last 10% of timesteps."
+        autofix_last || throw(ArgumentError("Cannot extract $last timesteps from a \
+                                            trajectory solution with only \
+                                            $(length(solution.t)) timesteps. \
+                                            Consider decreasing the `last` argument value \
+                                            and/or specifying it as a percentage instead \
+                                            (*e.g.* `'10%'`) or setting \
+                                            `autofix_last = true`."))
+
+        ten_perc = round(Int, n_timesteps * 10 / 100)
+        quiet || @warn "Cannot extract $last timesteps from a trajectory solution \
+        with only $(length(solution.t)) timesteps. \
+        We set 'last' to $(ten_perc), i.e. the last 10% of timesteps."
+        last = ten_perc
     end
 
-    out = solution[:, end-(last-1):end]
-
-    if !isnothing(idxs)
-        # Handle species names:
-        if any(eltype(idxs) .=== [String, Char])
-            sp = get_parameters(solution).network.species
-            if eltype(idxs) === Char
-                idxs = [idxs]
-            end
-            idxs = indexin(idxs, sp)
-        end
-        # Extract species indices:
-        out = out[idxs, :]
-    end
-    deepcopy(out)
+    last
 end
diff --git a/test/measures/test-functioning.jl b/test/measures/test-functioning.jl
index f612b5d35..16c57753f 100644
--- a/test/measures/test-functioning.jl
+++ b/test/measures/test-functioning.jl
@@ -1,3 +1,44 @@
+@testset "Extraction of living species" begin
+
+    foodweb = FoodWeb([0 0; 1 0])
+    params = ModelParameters(foodweb)
+    sol = simulates(params, [0.5, 0.5])
+
+
+    @test living_species(sol) == (species = ["s1", "s2"], idxs = [1, 2])
+    @test living_species(sol; idxs = 1) ==
+          living_species(sol; idxs = [1]) ==
+          living_species(sol; idxs = "s1") ==
+          living_species(sol; idxs = ["s1"]) ==
+          (species = ["s1"], idxs = [1])
+
+    @test living_species(sol; idxs = 2) ==
+          living_species(sol; idxs = [2]) ==
+          living_species(sol; idxs = "s2") ==
+          living_species(sol; idxs = ["s2"]) ==
+          (species = ["s2"], idxs = [2])
+
+    sol0 = simulates(params, [0, 0])
+
+    @test living_species(sol0).idxs == living_species([0 0; 0 0]) == Int64[]
+    @test living_species(sol0).species == String[]
+
+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)
+
+    @test_throws ArgumentError("`idxs` is not currently supported for \
+                               `trophic_structure()` and it may never be. \
+                               Please set `idxs = nothing`.") trophic_structure(
+        sol;
+        last = 10,
+        idxs = 2,
+    )
+end
+
 @testset "Producer growth rate" begin
 
     # Set up
@@ -98,7 +139,7 @@ end
           simpson(sim_one_sp; last = 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)) # Not defined for 0 species
+    @test isnan(simpson(sim_zero[:, end])) == isnan(simpson(sim_zero; last = 1))
 
     # Community evenness
     @test evenness(sim_two_sp[:, end]) == evenness(sim_two_sp; last = 1) ≈ 1.0 # Maximum equity of species biomass
diff --git a/test/measures/test-stability.jl b/test/measures/test-stability.jl
index 141e2d3c9..5e1eaa444 100644
--- a/test/measures/test-stability.jl
+++ b/test/measures/test-stability.jl
@@ -1,25 +1,36 @@
 @testset "temporal stability as CV" begin
 
+    mc = Matrix(undef, 0, 3)
+    mr = Matrix(undef, 3, 0)
 
     # Test avg Population stability
-    @test isnan(EcologicalNetworksDynamics.cv_sp([0 0; 0 0]).avg)
-    @test EcologicalNetworksDynamics.cv_sp([1 1; 1 1]).avg == 0
-    @test EcologicalNetworksDynamics.cv_sp([0 1; 1 0]).avg ==
-          EcologicalNetworksDynamics.cv_sp([0 0; 1 1]).avg ==
+    @test isnan(cv_sp([0 0; 0 0]).avg) ==
+          isnan(cv_sp([0 0; 0 0]).avg) ==
+          isnan(cv_sp(mc).avg) ==
+          isnan(cv_sp(mr).avg) ==
+          isnan(cv_sp(mc).sp) ==
+          isnan(cv_sp(mr).sp)
+    @test all(isnan.(cv_sp([0 0; 0 0]).sp))
+
+    @test cv_sp([1 1; 1 1]).avg == cv_sp([2 2; 2 2]).avg == 0
+    @test cv_sp([0 1; 1 0]).avg ==
+          cv_sp([0 0; 1 1]).avg ==
           std([0 1]; corrected = false) * 2 / (0.5 * 2)
 
     # Test synchrony
-    @test BEFWM2.synchrony([0 1; 1 0]) ≈ 0
-    @test BEFWM2.synchrony([0 0; 1 1]) ≈ 1
+    @test isnan(synchrony([0 0; 0 0])) == isnan(synchrony(mc)) == isnan(synchrony(mr))
+    @test synchrony([0 1; 1 0]) == synchrony([1 0; 0 1]) == 0.0
+    @test synchrony([0 0; 1 1]) == synchrony([0 0; 0 1]) == 1.0
 
     # Test CV decompostion in population statibility and synchrony
     mat = rand(10, 5)
     @test std(sum.(eachrow(mat)); corrected = false) / mean(sum.(eachrow(mat))) ≈
-          EcologicalNetworksDynamics.cv_sp(mat).avg *
-          sqrt(EcologicalNetworksDynamics.synchrony(mat))
-    @test EcologicalNetworksDynamics.temporal_cv(mat) ≈
-          EcologicalNetworksDynamics.cv_sp(mat).avg *
-          sqrt(EcologicalNetworksDynamics.synchrony(mat))
+          cv_sp(mat).avg * sqrt(synchrony(mat))
+    @test community_cv(mat) ≈ cv_sp(mat).avg * sqrt(synchrony(mat))
+
+    @test isnan(community_cv([0 0; 0 0])) ==
+          isnan(community_cv(mc)) ==
+          isnan(community_cv(mr))
 
     foodweb = FoodWeb([0 0; 1 0]; Z = 1) # Two producers and one co
     params = ModelParameters(
@@ -30,9 +41,37 @@
 
     cv_com = coefficient_of_variation(sol; last = 100)
     cv_one_sp = coefficient_of_variation(sol; idxs = [1], last = 100)
+    cv_one_sp2 = coefficient_of_variation(sol; idxs = [2], last = 100)
 
     # One species is fully synchronous with itself
-    @test cv_one_sp.synchrony ≈ 1 atol = 0.01
+    @test synchrony(sol; idxs = [1], last = 100) ==
+          cv_one_sp.synchrony ==
+          synchrony(sol; idxs = "s1", last = 100) ==
+          1.0
+    @test cv_one_sp2.synchrony ==
+          synchrony(sol; idxs = [2], last = 100) ==
+          synchrony(sol; idxs = "s2", last = 100)
+    @test cv_one_sp2.synchrony ≈ 1.0 atol = 10^-12
     # Strongly oscillating predator-prey are not so synchronous:
-    @test cv_com.synchrony <= 0.5
+    @test synchrony(sol; last = 100) == cv_com.synchrony <= 0.5
+
+    @test cv_one_sp.cv_sp[1] ==
+          cv_one_sp.avg_cv_sp ==
+          cv_one_sp.cv_com ==
+          cv_sp(sol; idxs = [1], last = 100).avg_cv_sp ==
+          cv_sp(sol; idxs = [1], last = 100).cv_sp[1] ==
+          cv_sp(sol; idxs = "s1", last = 100).avg_cv_sp ==
+          cv_sp(sol; idxs = "s1", last = 100).cv_sp[1]
+    @test cv_one_sp.cv_sp[1] ≈ 2.0 atol = 10^-1
+
+    # Test "dead" species removal with CV
+    one_sp = simulates(params, [0.25, 0.0]; tmax = 500, callback = nothing, t0 = 0)
+    cv_one_sp = coefficient_of_variation(one_sp; idxs = [1], last = 10)
+
+    @test synchrony(one_sp; last = 10) == cv_one_sp.synchrony == 1.0
+
+    cv_one_sp2 = coefficient_of_variation(one_sp; idxs = 2, last = 10)
+    @test all(map(x -> isnan(x), cv_one_sp2))
+
+
 end
diff --git a/test/measures/test-utils.jl b/test/measures/test-utils.jl
index d830a676d..9ca7f1e25 100644
--- a/test/measures/test-utils.jl
+++ b/test/measures/test-utils.jl
@@ -9,35 +9,95 @@
     @test round.(extract_last_timesteps(sim; last = "20.5%")) == [0.0 0.0; 1.0 1.0]
 
     # Test errors on sugar
-    @test_throws AssertionError("`last` Char should end by %") extract_last_timesteps(
+    @test_throws ArgumentError("`last` String should end by %.") extract_last_timesteps(
         sim,
-        last = "20",
+        last = "15%5",
     )
-    @test_throws AssertionError("`last` should be > 0 and <= 100") extract_last_timesteps(
+
+    @test_throws ArgumentError("Cannot extract -100.1% of the solution's timesteps: 0% < \
+                               `last` <= 100% must hold.") extract_last_timesteps(
+        sim,
+        last = "-100.1%",
+    )
+
+    @test_throws ArgumentError("Cannot extract 0.0% of the solution's timesteps: 0% < \
+                               `last` <= 100% must hold.") extract_last_timesteps(
         sim,
         last = "0%",
     )
-    @test_throws AssertionError("`last` should be > 0 and <= 100") extract_last_timesteps(
+
+    # Test last as integer
+    @test_throws ArgumentError(
+        "Cannot extract 100 timesteps from a trajectory solution \
+        with only 12 timesteps. Consider decreasing the `last` \
+        argument value and/or specifying it as a percentage instead \
+        (*e.g.* `'10%'`) or setting `autofix_last = true`.",
+    ) extract_last_timesteps(sim, last = 100)
+    @test_throws ArgumentError("Cannot extract 0 timesteps. `last` should be \
+                               an Integer strictly superior to 0.") extract_last_timesteps(
+        sim,
+        last = 0,
+    )
+    @test_throws ArgumentError("Cannot extract -10 timesteps. `last` should be \
+                               an Integer strictly superior to 0.") extract_last_timesteps(
         sim,
-        last = "100.1%",
+        last = -10,
     )
 
-    @test_throws AssertionError("`last` > length(solution.t),
-        please decrease `last` \
-        or consider to setting `autofix_last = true`.") extract_last_timesteps(
+    @test_throws ArgumentError("Cannot extract `last` from a decimal number. \
+                               Did you mean \"4.5%\"?") extract_last_timesteps(
         sim,
-        last = 100,
+        last = 4.5,
     )
-    @test_warn "'last' is lower or equal to the number of timesteps. \
-    'last' is then set to 1, \
-    i.e. the last 10% of timesteps." extract_last_timesteps(
+
+    @test_throws ArgumentError("Cannot extract timesteps with \
+                               Matrix{Any}(undef, 0, 0). \
+                               `last` should be an Integer or a String \
+                               representing a percentage.") extract_last_timesteps(
         sim,
-        last = 100,
-        autofix_last = true,
+        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)
     @test_nowarn extract_last_timesteps(sim, last = 100, autofix_last = true, quiet = true)
 
+    # Test species selection
+    @test extract_last_timesteps(sim; last = "10%", idxs = ["s1"]) ==
+          extract_last_timesteps(sim; last = "10%", idxs = "s1") ==
+          [0.0;;]
+
+    err_sp_msg = "Species [\"s3\"] are not found in the network. Any mispelling?"
+    @test_throws ArgumentError(err_sp_msg) extract_last_timesteps(
+        sim;
+        last = "10%",
+        idxs = "s3",
+    )
+    @test_throws ArgumentError(err_sp_msg) extract_last_timesteps(
+        sim;
+        last = "10%",
+        idxs = ["s3"],
+    )
+
+    err_id_msg = "Cannot extract idxs [4]. Each species index should  be \
+                               superior or equal to 1 and should be inferior or equal to \
+                               the number of species (2)."
+    @test_throws ArgumentError(err_id_msg) extract_last_timesteps(
+        sim;
+        last = "10%",
+        idxs = 4,
+    )
+
+    @test_throws ArgumentError(err_id_msg) extract_last_timesteps(
+        sim;
+        last = "10%",
+        idxs = [2, 4],
+    )
+
+
     @test size(extract_last_timesteps(sim; last = 1), 2) == 1
     @test size(extract_last_timesteps(sim; last = 10), 2) == 10
-    @test extract_last_timesteps(sim; last = 10) isa Matrix
+    @test extract_last_timesteps(sim; last = 10) isa AbstractMatrix
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index c5fdcdb6a..3d354ea0c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,7 +3,6 @@ using EcologicalNetworksDynamics
 using Test
 using SparseArrays
 using Random
-using EcologicalNetworks
 using JuliaFormatter
 using SyntaxTree
 using Logging #  TODO: remove once warnings are removed from `generate_dbdt`.
@@ -11,7 +10,7 @@ using Statistics
 
 
 # Set and print seed
-seed = sample(1:100000)
+seed = rand(1:100000)
 Random.seed!(seed)
 @info "Seed set to $seed."