WIP: Incorporate Iago's suggestions

use_cases/biomass-vs-mortality.jl

@@ -1,4 +1,3 @@
-# Import dependencies.
 using EcologicalNetworksDynamics
 using DifferentialEquations
 using DiffEqCallbacks
@@ -15,7 +14,7 @@ i_intra = 0.8 # intraspecific interference
 
 # Generate food webs.
 n_foodweb = 50
-foodweb_vector = [FoodWeb(nichemodel, S; C = C, tol = 0.01, Z = Z) for i in 1:n_foodweb]
+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
@@ -34,32 +33,33 @@ 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))
-        j != 1 || (A_ref = remove_extinct_species(fw.A, extinct_sp)) # update reference for d=0
-        class_matrix[i, j, :] = length.(collect(trophic_class(A_ref, extinct_sp = extinct_sp))) # save results
+        # 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
     end
 end
+println("All simulations done.") # conclusion
 
 # Process data.
 relative_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])
 end
-err = reshape(std(relative_class_variation; dims = 1), n_d, n_class) ./ sqrt(n_foodweb)
+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)
 
 # Plot.
-col = reshape(collect(palette(:viridis, n_class)), 1, n_class) # define our color palette
-plot(
-    d0_range,
-    reshape(mean(relative_class_variation; dims = 1), n_d, n_class);
-    label = "",
-    color = col,
-    lw = 2,
-)
+col = collect(palette(:viridis, n_class))' # define our color palette, (' takes the transpose)
+plot(d0_range, mean_variation; label = "", color = col, lw = 2)
 scatter!(
     d0_range,
-    reshape(mean(relative_class_variation; dims = 1), n_d, n_class);
+    mean_variation;
     labels = ["producers" "intermediate consumers" "top predators"],
-    yerr = err,
+    yerr = err_variation,
     markersize = 5,
     markerstrokewidth = 1.2,
     size = (500, 500),

use_cases/diversity-vs-nti-intensity.jl

@@ -1,4 +1,3 @@
-# Set up parallelization and import dependencies.
 using DiffEqCallbacks
 using DifferentialEquations
 using EcologicalNetworksDynamics
@@ -18,40 +17,34 @@ 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 = C, tol = 0.01, Z = Z) for i in 1:n_foodweb]
+foodweb_backbone = [FoodWeb(nichemodel, S; C, Z, tol = 0.01) for _ in 1:n_foodweb]
 
 # Non-trophic interactions.
-n_intensity = 5
+intensity_range_size = 5
+range_to(max) = LinRange(0, max, intensity_range_size)
 intensity_dict = Dict(
-    :interference => LinRange(0, 4, n_intensity),
-    :facilitation => LinRange(0, 4, n_intensity),
-    :refuge => LinRange(0, 5, n_intensity),
-    :competition => LinRange(0, 0.4, n_intensity),
+    :interference => range_to(4),
+    :facilitation => range_to(4),
+    :refuge => range_to(5),
+    :competition => range_to(0.4),
 )
-interaction_names = vcat(keys(intensity_dict)...)
+interaction_names = keys(intensity_dict)
 n_interaction = length(interaction_names)
-diversity = zeros(n_interaction, n_intensity, n_foodweb)
-"""
-Keyword arguments to pass the number of `L`inks and the `i`ntensity
-for the non-trophic interactions `nti_name`.
-"""
-function get_kwargs(nti_name, L, i) 
-    nti_name_str = String(nti_name)
-    Dict(Symbol("L_" * nti_name_str) => L, Symbol("I_" * nti_name_str) => i)
-end
+diversity = zeros(n_interaction, intensity_range_size, n_foodweb) # store diversity at equilibrium
 
 # 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>`, 
+# To run simulations in parallel one should execute the script as follow
+# `$ 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.
 Threads.@threads for i in 1:n_foodweb # parallelize if possible
     foodweb = foodweb_backbone[i]
     for (j, interaction) in enumerate(interaction_names)
         intensity_range = intensity_dict[interaction]
         for (k, intensity) in enumerate(intensity_range)
-            kwargs = get_kwargs(interaction, L_nti, intensity)
-            net = MultiplexNetwork(foodweb; kwargs...)
+            net = MultiplexNetwork(foodweb; [interaction => (L = L_nti, I = intensity)]...)
             F = ClassicResponse(net; c = i_intra)
             p = ModelParameters(net; functional_response = F)
             cb_set = CallbackSet(
@@ -61,34 +54,34 @@ 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 done.") # print progress
+        println("Interaction $interaction (foodweb $i) done.") # full info after every long simulation
     end
-    println("Food web $i done.") # print progress
 end
+println("All simulations done.") # conclusion
 
 # Compute variation in diversity between trophic only and
 # trophic with non-trophic interactions.
-diversity_variation = zeros(size(diversity))
-for i in 1:n_interaction, j in 1:n_intensity, k in 1:n_foodweb
-    diversity_variation[i, j, k] = (diversity[i, j, k] - diversity[i, 1, k]) / S
+variation = zeros(size(diversity))
+for i in 1:n_interaction, j in 1:intensity_range_size, k in 1:n_foodweb
+    variation[i, j, k] = (diversity[i, j, k] - diversity[i, 1, k]) / S
 end
-mean_diversity_variation = mean(diversity_variation; dims = 3)
-err_diversity_variation = std(diversity_variation; dims = 3) ./ sqrt(n_foodweb)
+mean_variation = mean(variation; dims = 3)
+err_variation = std(variation; dims = 3) ./ sqrt(n_foodweb)
 
 # Plot.
 plot_vec = []
 for (i, interaction) in enumerate(interaction_names)
     intensity_range = intensity_dict[interaction]
     p = scatter(
         intensity_range,
-        mean_diversity_variation[i, :];
-        yerr = err_diversity_variation[i, :],
+        mean_variation[i, :];
+        yerr = err_variation[i, :],
         color = :black,
         label = "",
     )
-    Interaction = uppercasefirst(String(interaction))
-    i0 = lowercase(first(String(interaction))) * "_0"
-    xlabel!(L"%$Interaction intensity, $%$i0$")
+    interaction = String(interaction)
+    interaction, i0 = uppercasefirst(interaction), first(interaction) * "_0"
+    xlabel!(L"%$interaction intensity, $%$i0$")
     ylabel!(L"Diversity variation, $\Delta S$")
     push!(plot_vec, p)
 end

use_cases/utils.jl

@@ -10,7 +10,7 @@ If that latter does not exist, it will be created.
 """
 function save_plot(fname_without_extension; verbose = true)
     checkdir("figures"; verbose = verbose)
-    fname_without_extension = "figures/" * fname_without_extension
+    fname_without_extension = joinpath("figures", fname_without_extension)
     format_list = ["png", "svg", "pdf"]
     for fmt in format_list
         fname_with_extension = fname_without_extension * "." * fmt
@@ -23,14 +23,14 @@ end
 """
     save_data(fname, data; verbose = false)
 
-Save `data` Julia object using `Serialization:serialize()` as `data/<fname>.jldata`.
-Data is automatically saved in the `data/` subdirectory and 
+Save `data` Julia object using `Serialization.serialize()` as `data/<fname>.jldata`.
+Data is automatically saved in the `data/` subdirectory and
 under the `.jldata` format.
 This choice is made to ensure consistency through the project.
 """
 function save_data(fname, data; verbose = true)
     checkdir("data"; verbose = verbose)
-    fname = "data/" * fname * ".jldata"
+    fname = joinpath("data", fname) * ".jldata"
     serialize(fname, data)
     !verbose || @info "$fname has been saved."
 end