Remove utils.jl, improve clarity, etc.

.gitignore

@@ -3,6 +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,149 @@ 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 in 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))
+    class = [prod, intermediate_cons, top_pred]
+    if !isempty(extinct_sp) # filter extinct species
+        class = [setdiff(cl, extinct_sp) for cl in class]
+    end
+    (producers = class[1], intermediate_consumers = class[2], top_predators = class[3])
+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/biomass-vs-mortality.jl

@@ -4,56 +4,80 @@ using DiffEqCallbacks
 using LaTeXStrings
 using Plots
 using Statistics
-include("utils.jl")
 
 # 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 
-
-# Generate food webs.
-n_foodweb = 50
-foodweb_vector = [FoodWeb(nichemodel, S; C, Z, tol = 0.01) for _ in 1:n_foodweb]
-
-n_d = 6
-d0_range = LinRange(0, 0.5, n_d) # mortality range
-n_class = 3 # 3 trophic classes: plants, intermediate consumers and top predators
+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
-    A_ref = nothing # where to store the reference trophic network
-    fw = foodweb_vector[i]
+    # 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)
-        F = ClassicResponse(fw; c = i_intra)
         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, rand(S); tmax = 10_000, callback = cb_set)
+        sol = simulate(
+            p,
+            B0;
+            tmax = 10_000,
+            callback = cb_set,
+            diff_code_data = (dBdt!, data),
+        )
         extinct_sp = keys(get_extinct_species(sol))
-        # Update reference for null mortality (d0 = 0).
-        j == 1 && (A_ref = remove_extinct_species(fw.A, extinct_sp))
-        # Record results.
-        class_matrix[i, j, :] =
-            length.(values(trophic_class(A_ref; extinct_sp = extinct_sp)))
-        println("d0 = $d0 (foodweb $i) done.") # full info after every long run
+        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.") # conclusion
+println("All simulations done.")
 
-# Process data.
-relative_class_variation = zeros(Float64, size(class_matrix))
+# 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
-    relative_class_variation[i, j, k] = 1 - (class_matrix[i, j, k] / class_matrix[i, 1, k])
+    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(relative_class_variation; dims = 1), n_d, n_class);
-err_variation =
-    reshape(std(relative_class_variation; dims = 1), n_d, n_class) ./ sqrt(n_foodweb)
+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.
-col = collect(palette(:viridis, n_class))' # define our color palette, (' takes the transpose)
+# 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,
@@ -67,13 +91,5 @@ scatter!(
 )
 xlabel!(L"Mortality rate, $d_0$")
 ylabel!(L"Relative loss of trophic class, $\gamma$")
-save_plot("loss-trophic-class-vs-mortality")
-
-# Save data.
-data_to_save = Dict(
-    :d0_range => d0_range,
-    :class_matrix => class_matrix,
-    :relative_class_variation => relative_class_variation,
-    :foodweb_vector => foodweb_vector,
-)
-save_data("trophic-sensitivity-vs-mortality", data_to_save)
+plot!(; size = (MathConstants.φ * 350, 350))
+# savefig("figures/loss-trophic-class-vs-mortality.pdf")