Refresh use_cases/ and introduce the first two.

- Mortality.
- Non-trophic interations intensity *vs.* diversity.

.gitignore

@@ -3,4 +3,5 @@ Manifest.toml
 LocalPreferences.toml
 docs/Manifest.toml
 docs/build/
-.DS_Store
+use_cases/figures/
+use_cases/data/

src/EcologicalNetworksDynamics.jl

@@ -145,13 +145,16 @@ export preys_of
 export producer_growth
 export ProducerCompetition
 export producers
+export remove_species
 export richness
 export set_temperature!
 export simulate
 export species_persistence
 export species_richness
 export TemperatureResponse
+export top_predators
 export total_biomass
+export trophic_classes
 export trophic_levels
 
 end

src/measures/structure.jl

@@ -1,6 +1,4 @@
-#=
-Quantifying food webs structural properties
-=#
+# Quantify food webs structural properties.
 
 """
 Number of species in the network.
@@ -182,7 +180,17 @@ function _ftl(A::AbstractMatrix)
 end
 
 """
+    trophic_levels(net::EcologicalNetwork)
+
 Trophic level of each species.
+
+```jldoctest
+julia> foodweb = FoodWeb([1 => 2])
+       trophic_levels(foodweb) == [2.0, 1.0]
+true
+```
+
+See also [`top_predators`](@ref) and [`trophic_classes`](@ref).
 """
 function trophic_levels(A::AbstractMatrix)
     A = Bool.(A)
@@ -192,8 +200,151 @@ function trophic_levels(A::AbstractMatrix)
     species_idx = parse.(Int64, [split(t, "s")[2] for t in trophiclvl_keys])
     trophiclvl_val[sortperm(species_idx)]
 end
+
 trophic_levels(net::EcologicalNetwork) = trophic_levels(get_trophic_adjacency(net))
 
+"""
+    top_predators(net::EcologicalNetwork)
+
+Top predator indexes (species eaten by nobody) of the given `net`work.
+
+```jldoctest
+julia> foodweb = FoodWeb([1 => 2, 2 => 3])
+       top_predators(foodweb) == [1]
+true
+```
+
+See also [`trophic_levels`](@ref) and [`trophic_classes`](@ref).
+"""
+function top_predators(net::EcologicalNetwork)
+    A = get_trophic_adjacency(net)
+    top_predators(A)
+end
+
+function top_predators(A)
+    S = richness(A)
+    [i for i in 1:S if all(A[:, i] .== 0)]
+end
+
+"""
+    trophic_classes(foodweb::FoodWeb)
+
+Categorize each species into three classes:
+either `producers`, `intermediate_consumers` or `top_predators`.
+
+```jldoctest
+julia> foodweb = FoodWeb([1 => 2, 2 => 3])
+       trophic_classes(foodweb)
+(producers = [3], intermediate_consumers = [2], top_predators = [1])
+```
+
+To access the name of the trophic classes
+you can simply call `trophic_classes` without any argument.
+
+```jldoctest
+julia> trophic_classes()
+(:producers, :intermediate_consumers, :top_predators)
+```
+
+See also [`trophic_levels`](@ref) and [`top_predators`](@ref).
+"""
+function trophic_classes(A; extinct_sp = Set())
+    S = richness(A)
+    prod = producers(A)
+    top_pred = top_predators(A)
+    intermediate_cons = setdiff(1:S, union(prod, top_pred))
+    # Order here must match the semantics of `trophic_classes()` order.
+    class = [prod, intermediate_cons, top_pred]
+    if !isempty(extinct_sp) # filter extinct species
+        class = [setdiff(cl, extinct_sp) for cl in class]
+    end
+    names = trophic_classes()
+    NamedTuple{names}(class)
+end
+
+trophic_classes() = (:producers, :intermediate_consumers, :top_predators)
+
+"""
+    remove_species(net::EcologicalNetwork, species_to_remove)
+
+Remove a list of species from a `FoodWeb`, a `MultiplexNetwork` or an adjacency matrix.
+
+# Examples
+
+Adjacency matrix.
+
+```jldoctest
+julia> A = [0 0 0; 1 0 0; 1 0 0]
+       remove_species(A, [1]) == [0 0; 0 0]
+true
+
+julia> remove_species(A, [2]) == [0 0; 1 0]
+true
+```
+
+`FoodWeb`.
+
+```jldoctest
+julia> foodweb = FoodWeb([0 0 0; 1 0 0; 1 0 0]; M = [1, 10, 20])
+       new_foodweb = remove_species(foodweb, [2])
+       new_foodweb.A == [0 0; 1 0]
+true
+
+julia> new_foodweb.M == [1, 20]
+true
+```
+
+`MultiplexNetwork`.
+
+```jldoctest
+julia> foodweb = FoodWeb([0 0 0; 1 0 0; 1 0 0]; M = [1, 10, 20])
+       net = MultiplexNetwork(foodweb; A_facilitation = [0 0 0; 1 0 0; 0 0 0])
+       new_net = remove_species(net, [2])
+       new_net.layers[:trophic].A == [0 0; 1 0]
+true
+
+julia> new_net.layers[:facilitation].A == [0 0; 0 0]
+true
+
+julia> new_net.M == [1, 20]
+true
+```
+"""
+function remove_species(foodweb::FoodWeb, species_to_remove)
+    S = richness(foodweb)
+    A_new = remove_species(foodweb.A, species_to_remove)
+    species_to_keep = setdiff(1:S, species_to_remove)
+    M_new = foodweb.M[species_to_keep]
+    metabolic_class_new = foodweb.metabolic_class[species_to_keep]
+    species_new = foodweb.species[species_to_keep]
+    method = foodweb.method
+    FoodWeb(A_new, species_new, M_new, metabolic_class_new, method)
+end
+
+function remove_species(net::MultiplexNetwork, species_to_remove)
+    trophic = convert(FoodWeb, net)
+    new_trophic = remove_species(trophic, species_to_remove)
+    multi_net = MultiplexNetwork(new_trophic) # no non-trophic interactions yet
+    for (name, layer) in net.layers
+        A_nti_new = remove_species(layer.A, species_to_remove)
+        multi_net.layers[name].A = A_nti_new
+        multi_net.layers[name].intensity = layer.intensity
+        multi_net.layers[name].f = layer.f
+    end
+    multi_net
+end
+
+function remove_species(A::AbstractMatrix, species_to_remove)
+    S = richness(A)
+    species_to_keep = setdiff(1:S, species_to_remove)
+    n_kept = length(species_to_keep)
+    A_simplified = zeros(Integer, n_kept, n_kept)
+    for (i_idx, i) in enumerate(species_to_keep), (j_idx, j) in enumerate(species_to_keep)
+        A_simplified[i_idx, j_idx] = A[i, j]
+    end
+    A_simplified
+end
+
 function massratio(obj::Union{ModelParameters,FoodWeb})
 
     if isa(obj, ModelParameters)

use_cases/Project.toml

@@ -1,7 +1,8 @@
 [deps]
-DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-EcologicalNetworks = "f03a62fe-f8ab-5b77-a061-bb599b765229"
-Mangal = "b8b640a6-63d9-51e6-b784-5033db27bef2"
+DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+EcologicalNetworksDynamics = "2fd9189a-c387-4076-88b7-22b33b5a4388"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781"
-PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
+Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

use_cases/biomass-vs-mortality.jl

@@ -0,0 +1,94 @@
+using EcologicalNetworksDynamics
+using DifferentialEquations
+using DiffEqCallbacks
+using LaTeXStrings
+using Plots
+using Statistics
+
+# Define global community parameters.
+S = 50 # species richness
+C = 0.06 # trophic connectance
+Z = 50 # predator-prey bodymass ratio
+i_intra = 0.8 # intraspecific interference
+n_foodweb = 50 # number of replicates of trophic networks
+n_d = 6 # number of tested mortality values
+d0_range = LinRange(0, 1, n_d) # mortality range
+n_class = length(trophic_classes()) # plants, intermediate consumers and top predators
+class_matrix = zeros(n_foodweb, n_d, n_class) # store the diversity in each trophic class
+
+# Main simulation loop.
+Threads.@threads for i in 1:n_foodweb # parallelize computation when possible
+    # We first define the food web structure as well as the trophic parameters
+    # that are kept constant. We check that there is no cycle so that the trophic levels
+    # and thus the trophic classes are well defined.
+    fw = FoodWeb(nichemodel, S; C, Z, check_cycle = true)
+    F = ClassicResponse(fw; c = i_intra) # set intraspecific predator interference
+    B0 = rand(S) # initial biomass associated with the food web
+    classes = nothing
+    # Then we loop on the mortality intensity that we increase at each iteration.
+    for (j, d0) in enumerate(d0_range)
+        br = BioRates(fw; d = d0 .* allometric_rate(fw, DefaultMortalityParams()))
+        p = ModelParameters(fw; functional_response = F, biorates = br)
+        # Prepare simulation boost.
+        dBdt_expr, data = generate_dbdt(p, :compact)
+        dBdt! = eval(dBdt_expr)
+        # For this specific set-up we need to define callbacks manually
+        # to ensure that the simulation stops when a fixed point is reached,
+        # and not before.
+        # To do so, we have lowered the `abstol` and `reltol` arguments of
+        # `TerminateSteadyState`.
+        # Moreover, we have set the `verbose` argument of the `ExtinctionCallback`
+        # to remove printed information about species extinctions.
+        cb_set =
+            CallbackSet(ExtinctionCallback(1e-5, false), TerminateSteadyState(1e-8, 1e-6))
+        sol = simulate(
+            p,
+            B0;
+            tmax = 10_000,
+            callback = cb_set,
+            diff_code_data = (dBdt!, data),
+        )
+        extinct_sp = keys(get_extinct_species(sol))
+        surviving_sp = setdiff(1:richness(fw), extinct_sp)
+        # Species classes are given by the the dynamics with zero mortality rate.
+        d0 == 0 && (classes = trophic_classes(remove_species(fw, extinct_sp)))
+        # Then record the surviving species in each class
+        # after running the dynamic with an increased mortality.
+        for k in 1:n_class
+            class_matrix[i, j, k] = length(intersect(classes[k], surviving_sp))
+        end
+        println("d0 = $d0 (foodweb $i) done.")
+    end
+end
+println("All simulations done.")
+
+# We compute the relative loss in each trophic class
+# for the different mortality intensities
+# compared to the classes for a zero mortality rate.
+class_variation = zeros(Float64, size(class_matrix))
+for i in 1:n_foodweb, j in 1:n_d, k in 1:n_class
+    if class_matrix[i, 1, k] == 0
+        class_variation[i, j, k] = 0
+    else
+        class_variation[i, j, k] = 1 - (class_matrix[i, j, k] / class_matrix[i, 1, k])
+    end
+end
+mean_variation = reshape(mean(class_variation; dims = 1), n_d, n_class);
+err_variation = reshape(std(class_variation; dims = 1), n_d, n_class) ./ sqrt(n_foodweb)
+
+# Plot relative loss in each trophic class versus the mortality rate intensity.
+col = collect(palette(:viridis, n_class))'
+plot(d0_range, mean_variation; label = "", color = col, lw = 2)
+scatter!(
+    d0_range,
+    mean_variation;
+    labels = ["producers" "intermediate consumers" "top predators"],
+    yerr = err_variation,
+    markersize = 5,
+    markerstrokewidth = 1.2,
+    size = (500, 500),
+    color = col,
+)
+xlabel!(L"Mortality rate, $d_0$")
+ylabel!(L"Relative loss of trophic class, $\gamma$")
+plot!(; size = (MathConstants.φ * 350, 350))