Remove utils.jl, improve clarity, etc.

src/EcologicalNetworksDynamics.jl

@@ -121,13 +121,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,139 @@ 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])
+```
+
+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
+
+"""
+    remove_species(net::EcologicalNetwork, species_to_remove)
+
+Remove a list of species from a `FoodWeb`, a `MultiplexNetwork` or an adjacency matrix.
+
+# Examples
+
+Adjancecy 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,27 +4,26 @@ 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 
+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
+d0_range = LinRange(0, 1, n_d) # mortality range
 n_class = 3 # 3 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]
+    classes = nothing
     for (j, d0) in enumerate(d0_range)
         F = ClassicResponse(fw; c = i_intra)
         br = BioRates(fw; d = d0 .* allometric_rate(fw, DefaultMortalityParams()))
@@ -33,11 +32,14 @@ Threads.@threads for i in 1:n_foodweb # parallelize computation when possible
             CallbackSet(ExtinctionCallback(1e-5, false), TerminateSteadyState(1e-8, 1e-6))
         sol = simulate(p, rand(S); tmax = 10_000, callback = cb_set)
         extinct_sp = keys(get_extinct_species(sol))
+        surviving_sp = setdiff(1:richness(fw), extinct_sp)
         # Update reference for null mortality (d0 = 0).
-        j == 1 && (A_ref = remove_extinct_species(fw.A, extinct_sp))
+        d0 == 0 && (classes = trophic_classes(remove_species(fw, extinct_sp)))
+        println(surviving_sp)
         # Record results.
-        class_matrix[i, j, :] =
-            length.(values(trophic_class(A_ref; extinct_sp = extinct_sp)))
+        for k in 1:n_class
+            class_matrix[i, j, k] = length(intersect(classes[k], surviving_sp))
+        end
         println("d0 = $d0 (foodweb $i) done.") # full info after every long run
     end
 end
@@ -67,13 +69,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")

use_cases/diversity-vs-nti-intensity.jl

@@ -4,7 +4,6 @@ using EcologicalNetworksDynamics
 using LaTeXStrings
 using Plots
 using Statistics
-include("utils.jl") # utility functions
 
 # Define global community parameters.
 S = 50 # species richness
@@ -17,25 +16,25 @@ L_nti += isodd(L_nti) # ensure that number of link is even for symmetric interac
 
 # Generate trophic backbones.
 n_foodweb = 10
-foodweb_backbone = [FoodWeb(nichemodel, S; C, Z, tol = 0.01) for _ in 1:n_foodweb]
+foodweb_backbone = [FoodWeb(nichemodel, S; C, Z) for _ in 1:n_foodweb]
 
-# Non-trophic interactions.
+# Set up non-trophic interactions.
 intensity_range_size = 5
 range_to(max) = LinRange(0, max, intensity_range_size)
-intensity_dict = Dict(
+intensity_dict = Dict( # intensity range for each non-trophic interaction
     :interference => range_to(4),
     :facilitation => range_to(4),
     :refuge => range_to(5),
     :competition => range_to(0.4),
 )
 interaction_names = keys(intensity_dict)
 n_interaction = length(interaction_names)
-diversity = zeros(n_interaction, intensity_range_size, n_foodweb) # store diversity at equilibrium
+diversity = zeros(n_interaction, intensity_range_size, n_foodweb) # store final diversity
 
 # Main simulation loop.
 # We parallelize the outer loop (when possible) with the `@threads` macro.
 # To run simulations in parallel one should execute the script as follow
-# `$ julia --threads <number_of_threads> <script_name.jl>`, 
+# `$ julia --threads <number_of_threads> <script_name.jl>`,
 # or open the Julia REPL with
 # `$ julia --threads <number_of_threads>`
 # where <number_of_threads> should be replace with the value of your choice.
@@ -54,7 +53,7 @@ Threads.@threads for i in 1:n_foodweb # parallelize if possible
             sol = simulate(p, rand(S); tmax = 10_000, callback = cb_set, verbose = false)
             diversity[j, k, i] = species_richness(sol[end])
         end
-        println("Interaction $interaction (foodweb $i) done.") # full info after every long simulation
+        println("Interaction $interaction (foodweb $i) done.") # info at each iteration
     end
 end
 println("All simulations done.") # conclusion
@@ -86,4 +85,3 @@ for (i, interaction) in enumerate(interaction_names)
     push!(plot_vec, p)
 end
 plot(plot_vec...)
-save_plot("nti-intensity-vs-diversity")