Change default threshold to 0

src/measures/functioning.jl

@@ -4,7 +4,7 @@ Adapted from BioenergeticFoodWeb.jl
 =#
 
 """
-    richness(solution; threshold = eps(), kwargs...)
+    richness(solution; threshold = 0, kwargs...)
 
 Returns the average number of species with a biomass larger than `threshold`
 over the `last` timesteps. `kwargs...` are optional arguments passed to
@@ -13,7 +13,9 @@ over the `last` timesteps. `kwargs...` are optional arguments passed to
 # Arguments:
 
   - `solution`: output of `simulate()` or `solve()`
-  - `threshold`: biomass threshold below which a species is considered extinct.
+  - `threshold`: biomass threshold below which a species is considered extinct. Set to 0 by
+  debault and it is recommended to let as this. It is recommended to change the threshold
+  using [`ExtinctionCallback`](@ref) in [`simulate`](@ref).
 
 # Examples
 
@@ -38,14 +40,14 @@ julia> even = evenness(sol);
 0.996
 ```
 """
-function richness(solution; threshold = eps(), kwargs...)
+function richness(solution; threshold = 0, kwargs...)
     measure_on = extract_last_timesteps(solution; kwargs...)
     rich = richness.(eachcol(measure_on), threshold = threshold)
     mean(rich)
 end
 
 """
-    richness(n::AbstractVector; threshold = eps())
+    richness(n::AbstractVector; threshold = 0)
 
 When applied to a vector of biomass, returns the number of biomass above `threshold`
 
@@ -59,7 +61,7 @@ julia> richness([1, 1])
 2
 ```
 """
-richness(n::AbstractVector; threshold = eps()) = sum(n .> threshold)
+richness(n::AbstractVector; threshold = 0) = sum(n .> threshold)
 
 """
     species_persistence(solution; threshold, kwargs...)
@@ -71,7 +73,7 @@ over the `last` timesteps.
 [`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
+`species_persistence(n::Vector; threshold = 0)`, it returns the proportion of
 species which biomass is above `threshold`.
 
 # Examples
@@ -103,7 +105,7 @@ function species_persistence(solution; kwargs...)
     m = richness(get_parameters(solution).network)
     r / m
 end
-species_persistence(n::AbstractVector; threshold = eps()) =
+species_persistence(n::AbstractVector; threshold = 0) =
     richness(n; threshold = threshold) / length(n)
 
 """
@@ -159,21 +161,21 @@ over the `last` timesteps.
 `kwargs...` arguments are forwarded to [`extract_last_timesteps`](@ref). See
 [`extract_last_timesteps`](@ref) for the argument details.
 
-Can also handle a vector, e.g. shannon_diversity(n::AbstractVector; threshold = eps())
+Can also handle a vector, e.g. shannon_diversity(n::AbstractVector; threshold = 0)
 
 # Reference
 
 https://en.wikipedia.org/wiki/Diversity_index#Shannon_index
 """
-function shannon_diversity(solution; threshold = eps(), kwargs...)
+function shannon_diversity(solution; threshold = 0, kwargs...)
     measure_on = extract_last_timesteps(solution; kwargs...)
 
     shan = shannon_diversity.(eachcol(measure_on), threshold = threshold)
     mean(shan)
 
 end
 
-function shannon_diversity(n::AbstractVector; threshold = eps())
+function shannon_diversity(n::AbstractVector; threshold = 0)
     x = filter(k -> k > threshold, n)
 
     if length(x) >= 1
@@ -196,20 +198,20 @@ over the `last` timesteps.
 `kwargs...` arguments are forwarded to [`extract_last_timesteps`](@ref). See
 [`extract_last_timesteps`](@ref) for the argument details.
 
-Can also handle a vector, e.g. simpson(n::AbstractVector; threshold = eps())
+Can also handle a vector, e.g. simpson(n::AbstractVector; threshold = 0)
 
 # Reference
 
 https://en.wikipedia.org/wiki/Diversity_index#Simpson_index
 """
-function simpson(solution; threshold = eps(), kwargs...)
+function simpson(solution; threshold = 0, kwargs...)
     measure_on = extract_last_timesteps(solution; kwargs...)
 
     simp = simpson.(eachcol(measure_on), threshold = threshold)
     mean(simp)
 end
 
-function simpson(n::AbstractVector; threshold = eps())
+function simpson(n::AbstractVector; threshold = 0)
     x = filter((k) -> k > threshold, n)
 
     if length(x) >= 1
@@ -230,20 +232,20 @@ Computes the average Pielou evenness, over the `last` timesteps.
 `kwargs...` arguments are forwarded to [`extract_last_timesteps`](@ref). See
 [`extract_last_timesteps`](@ref) for the argument details.
 
-Can also handle a vector, e.g. `evenness(n::AbstractVector; threshold = eps())`
+Can also handle a vector, e.g. `evenness(n::AbstractVector; threshold = 0)`
 
 # Reference
 
 https://en.wikipedia.org/wiki/Species_evenness
 """
-function evenness(solution; threshold = eps(), kwargs...)
+function evenness(solution; threshold = 0, kwargs...)
     measure_on = extract_last_timesteps(solution; kwargs...)
 
     piel = evenness.(eachcol(measure_on), threshold = threshold)
     mean(piel)
 end
 
-function evenness(n::AbstractVector; threshold = eps())
+function evenness(n::AbstractVector; threshold = 0)
     x = filter((k) -> k > threshold, n)
     if length(x) > 0
         even = shannon_diversity(x) / log(length(x))
@@ -370,7 +372,7 @@ julia> foodweb = FoodWeb([0 0 0; 0 1 0; 1 1 0]; quiet = true);
 -2
 ```
 """
-function trophic_structure(solution; threshold = eps(), idxs = nothing, kwargs...)
+function trophic_structure(solution; threshold = 0, idxs = nothing, kwargs...)
 
     isnothing(idxs) || throw(ArgumentError("`idxs` is not currently supported for \
                                            `trophic_structure()` and it may never be. \
@@ -424,7 +426,7 @@ julia> B0 = [0, 0.5, 0.5];
 (species = ["s2", "s3"], idxs = [2, 3])
 ```
 """
-function living_species(solution; threshold = eps(), idxs = nothing, kwargs...)
+function living_species(solution; threshold = 0, idxs = nothing, kwargs...)
 
     measure_on = extract_last_timesteps(solution; idxs = idxs, kwargs...)
 
@@ -439,8 +441,8 @@ function living_species(solution; threshold = eps(), idxs = nothing, kwargs...)
     (species = sp_alive_name, idxs = alive_sp)
 end
 
-living_species(mat::Matrix; threshold = eps()) =
-    findall(x -> x >= threshold, biomass(mat).sp)
+living_species(mat::Matrix; threshold = 0) =
+    findall(x -> x > threshold, biomass(mat).sp)
 
 """
     min_max(solution; kwargs...)

src/measures/stability.jl

@@ -9,7 +9,7 @@ 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...)
+function cv_sp(solution; threshold = 0, kwargs...)
 
     measure_on = extract_last_timesteps(solution; kwargs...)
     # Fetch species that are alive, whose avg biomass is > threshold
@@ -55,7 +55,7 @@ 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(solution; threshold = eps(), kwargs...)
+function synchrony(solution; threshold = 0, kwargs...)
 
     measure_on = extract_last_timesteps(solution; kwargs...)
     # Fetch species that are alive, whose avg biomass is > threshold
@@ -91,7 +91,7 @@ Compute the temporal Coefficient of Variation of community biomass.
 
 See [`coefficient_of_variation`](@ref) for the argument details.
 """
-function community_cv(solution; threshold = eps(), kwargs...)
+function community_cv(solution; threshold = 0, kwargs...)
 
     measure_on = extract_last_timesteps(solution; kwargs...)
     # Fetch species that are alive, whose avg biomass is > threshold
@@ -150,7 +150,7 @@ julia> B0 = [0, 0.5, 0.5]; # Two producers
 ((:cv_com, :avg_cv_sp, :synchrony), (0.14, 0.14, 1.0))
 ```
 """
-function coefficient_of_variation(solution; threshold = eps(), kwargs...)
+function coefficient_of_variation(solution; threshold = 0, kwargs...)
 
     measure_on = extract_last_timesteps(solution; kwargs...)
     # Fetch species that are alive, whose avg biomass is > threshold

src/measures/utils.jl

@@ -45,11 +45,15 @@ function extract_last_timesteps(solution; idxs = nothing, kwargs...)
 
     last = process_last_timesteps(solution; kwargs...)
 
+
     out = solution[:, end-(last-1):end]
 
     idxs = process_idxs(solution, idxs;)
     # Extract species indices:
     out = out[idxs, :]
+
+    quiet || check_last_extinction(solution; idxs = idxs, last = last)
+
     deepcopy(out)
 end
 
@@ -136,3 +140,71 @@ function process_last_timesteps(solution; last = "10%", quiet = false)
                                               (e.g. `\"10%\"`)."))
     last
 end
+
+function check_last_extinction(solution; idxs = nothing, last)
+    ext = get_extinction_timesteps(solution; idxs)
+    ext_t = ext.extinction_timestep
+
+    extinct = ext_t .!== nothing
+    if any(extinct)
+        check_last = findall(x -> x > last, ext_t[extinct])
+
+        isempty(check_last) || @warn "With `last` = $last, a table has been extracted with \
+        the extinct species $(ext.species[extinct][check_last]), \
+        which went extinct at timesteps = $(ext_t[extinct][check_last]). Set `last` >= \
+        $(maximum(ext_t[extinct])) to get rid of them."
+    end
+end
+
+function get_extinction_timesteps(solution; idxs = nothing)
+    idxs = process_idxs(solution, idxs;)
+    sp = get_parameters(solution).network.species[idxs]
+
+    ext_t = findfirst.(isequal(0), eachrow(solution[idxs,:]))
+    extinct = ext_t .!== nothing
+
+    (
+     species = sp[extinct],
+     idxs = idxs[extinct],
+     extinction_timestep =  something.(ext_t[extinct])
+    )
+end
+
+function get_extinction_timesteps(m::AbstractVector; threshold = 0)
+     findfirst(x -> x <= threshold, m)
+end
+
+
+function get_extinction_timesteps(solution; idxs = nothing)
+    idxs = process_idxs(solution, idxs;)
+    sp = get_parameters(solution).network.species[idxs]
+
+    ext_t = findfirst.(isequal(0), eachrow(solution[idxs,:]))
+    extinct = ext_t .!== nothing
+
+    (
+     species = sp[extinct],
+     idxs = idxs[extinct],
+     extinction_timestep =  ext_t[extinct]
+    )
+end
+
+function get_alive_species(solution; idxs = nothing, threshold = 0)
+    idxs = process_idxs(solution, idxs;)
+    sp = get_parameters(solution).network.species[idxs]
+
+    alive_idxs = get_alive_species(solution[idxs, end]; threshold = threshold)
+
+    (
+     species = sp[alive_idxs],
+     idxs = idxs[alive_idxs]
+    )
+end
+
+function get_extinct_species(m::AbstractVector; threshold = 0)
+    findall(x -> x <= threshold, m)
+end
+
+function get_alive_species(m::AbstractVector; threshold = 0)
+    findall(x -> x > threshold, m)
+end