Integrate Iago's suggestions

use_cases/diversity-vs-nti-intensity.jl

@@ -23,11 +23,9 @@ n_foodweb = 50 # Replicates of trophic backbones.
 # 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.
 tmax = 10_000
-callback = CallbackSet(ExtinctionCallback(1e-5, false), TerminateSteadyState(1e-8, 1e-6))
-verbose = false # Do not show '@info' message when species go extinct.
+verbose = false # Do not show '@info' messages during simulation run.
+callback = CallbackSet(ExtinctionCallback(1e-5, verbose), TerminateSteadyState(1e-8, 1e-6))
 
 # Set up non-trophic interactions.
 intensity_range_size = 5
@@ -42,59 +40,56 @@ interaction_names = keys(intensity_values_dict)
 n_interaction = length(interaction_names)
 
 # Main simulation loop.
-# Each thread writes in its own DataFrame. We merge them at the end of the loop.
-df_to_merge = DataFrame[]
-Threads.@threads for i in 1:n_foodweb # We parallelize computation if possible.
-    df_tmp = DataFrame(;
+# Each thread handle a single food web and writes in its own DataFrame.
+# Merge the DataFrames at the end of the loop.
+dfs = [DataFrame() for _ in 1:n_foodweb] # Fill the vector with empty DataFrames.
+Threads.@threads for i in 1:n_foodweb # Parallelize computation if possible.
+    df_thread = DataFrame(;
         foodweb_id = Int64[],
         interaction = Symbol[],
         intensity = Float64[],
         diversity = Int64[],
     )
-    # First we create a trophic backbone.
+    # First, create a trophic backbone.
     foodweb = FoodWeb(nichemodel, S; C, Z)
-    # Secondly we loop on the different non-trophic interactions we want to study.
+    # Secondly, loop on the different non-trophic interactions to study.
     for interaction in interaction_names
         intensity_values = intensity_values_dict[interaction]
-        # Third, we increase the non-trophic interaction intensity.
+        # Third, increase the non-trophic interaction intensity.
         for intensity in intensity_values
-            # We add the given non-trophic interaction to the trophic backbone
+            # Add the given non-trophic interaction to the trophic backbone
             # with the given intensity.
             net = MultiplexNetwork(foodweb; [interaction => (L = L_nti, I = intensity)]...)
             functional_response = ClassicResponse(net; c = i_intra)
             params = ModelParameters(net; functional_response)
-            solution = simulate(params, rand(S); tmax, callback, verbose)
-            # We save the final diversity at equilibrium.
-            push!(df_tmp, [i, interaction, intensity, richness(solution[end])])
+            solution = simulate(params, rand(S); tmax, callback)
+            # Save the final diversity at equilibrium.
+            push!(df_thread, [i, interaction, intensity, richness(solution[end])])
         end
         @info "Interaction $interaction (foodweb $i) done."
     end
-    push!(df_to_merge, df_tmp)
+    dfs[i] = df_thread
 end
-df = reduce(vcat, df_to_merge)
 @info "All simulations done."
+df = reduce(vcat, dfs)
 
-# First, for each trophic backbone we compute the relative diversity difference
-# compared to the reference that is given by the trophic network that does not
-# contain non-trophic interactions, i.e. intensity = 0.
+# Compute the relative diversity different for each network
+# compared to the reference network that does not contain non-trophic interactions,
+# i.e. non-trophic intensity is null.
 groups = groupby(df, [:foodweb_id, :interaction])
-relative_difference(x, x_ref) = (x - x_ref) / x_ref
 for sub_df in groups
-    reference_row = findfirst(x -> x == 0, sub_df.intensity) # Trophic interactions only.
+    reference_row = findfirst(==(0), sub_df.intensity) # Trophic interactions only.
     S_ref = sub_df[reference_row, :diversity]
-    transform!(
-        sub_df,
-        :diversity => ByRow(S -> relative_difference(S, S_ref)) => :delta_diversity,
-    )
+    relative_difference(S) = (S - S_ref) / S_ref
+    transform!(sub_df, :diversity => ByRow(relative_difference) => :delta_diversity)
 end
-df = combine(groups, :) # Secondly, we recombine the DataFrames.
-# Third, for each interaction and intensity combination we compute the mean
-# and the confidence interval of diversity variation.
-ic95(x) = 1.96 * std(x) / sqrt(length(x)) # Confidence interval at 95%.
+df = combine(groups, :) # Recombine DataFrames.
+# Compute the mean and the confidence interval of diversity variation for each interaction.
+ci95(x) = 1.96 * std(x) / sqrt(length(x)) # Confidence interval at 95%.
 df_processed = combine(
     groupby(df, [:interaction, :intensity]),
     :delta_diversity => mean,
-    :delta_diversity => ic95,
+    :delta_diversity => ci95,
 )
 
 # Plot diversity variation versus non-trophic interaction intensity.
@@ -115,9 +110,9 @@ for (idx, interaction) in zip(indices, interaction_names)
     )
     x = df_interaction.intensity
     y = df_interaction.delta_diversity_mean
-    ic = df_interaction.delta_diversity_ic95
+    ci = df_interaction.delta_diversity_ci95
     scatterlines!(x, y)
-    errorbars!(x, y, ic; whiskerwidth = 5)
+    errorbars!(x, y, ci; whiskerwidth = 5)
 end
 
 # Write x and y labels.

use_cases/paradox-enrichment.jl

@@ -0,0 +1,60 @@
+import AlgebraOfGraphics: set_aog_theme!
+
+using CairoMakie
+using DataFrames
+using EcologicalNetworksDynamics
+
+"""
+Compute biomass extrema for each species during the `last` time steps.
+"""
+function biomass_extrema(solution, last)
+    trajectories = extract_last_timesteps(solution; last, quiet = true)
+    S = size(trajectories, 1) # Row = species, column = time steps.
+    [(min = minimum(trajectories[i, :]), max = maximum(trajectories[i, :])) for i in 1:S]
+end
+
+foodweb = FoodWeb([0 0; 1 0]); # 2 eats 1.
+functional_response = ClassicResponse(foodweb; aᵣ = 1, hₜ = 1, h = 1);
+K_values = LinRange(1, 10, 50)
+tmax = 1_000 # Simulation length.
+verbose = false # Do not show '@info' messages during the simulation.
+df = DataFrame(;
+    K = Float64[],
+    B_resource_min = Float64[],
+    B_resource_max = Float64[],
+    B_consumer_min = Float64[],
+    B_consumer_max = Float64[],
+)
+
+# Run simulations: compute equilibirum biomass for each carrying capacity.
+for K in K_values
+    environment = Environment(foodweb; K)
+    params = ModelParameters(foodweb; functional_response, environment)
+    B0 = rand(2) # Inital biomass.
+    solution = simulate(params, B0; tmax, verbose)
+    extrema = biomass_extrema(solution, "10%")
+    push!(df, [K, extrema[1].min, extrema[1].max, extrema[2].min, extrema[2].max])
+end
+
+# Plot the orbit diagram with Makie.
+set_aog_theme!() # AlgebraOfGraphics theme.
+c_r = :green # Resource color.
+c_c = :purple # Consumer color.
+c_v = :grey # Vertical lines color.
+fig = Figure()
+ax = Axis(fig[2, 1]; xlabel = "Carrying capacity, K", ylabel = "Equilibrium biomass")
+resource_line = scatterlines!(df.K, df.B_resource_min; color = c_r, markercolor = c_r)
+scatterlines!(df.K, df.B_resource_max; color = c_r, markercolor = c_r)
+consumer_line = scatterlines!(df.K, df.B_consumer_min; color = c_c, markercolor = c_c)
+scatterlines!(df.K, df.B_consumer_max; color = c_c, markercolor = c_c)
+K0 = 2.3
+v_line1 = vlines!(ax, [K0]; color = c_v)
+v_line2 = vlines!(ax, [1 + 2 * K0]; color = c_v, linestyle = :dashdot)
+Legend(
+    fig[1, 1],
+    [resource_line, consumer_line, v_line1, v_line2],
+    ["resource", "consumer", "K₀", "1+2K₀"];
+    orientation = :horizontal,
+    tellheight = true, # Adjust the height of the legend sub-figure.
+    tellwidth = false, # Do not adjust width of the orbit diagram.
+)

use_cases/persistence-vs-producer-competition.jl

@@ -1,9 +1,10 @@
-using AlgebraOfGraphics
+import AlgebraOfGraphics: set_aog_theme!
+import ColorSchemes: get, viridis
+import Statistics: mean, std
+
 using CairoMakie
-using ColorSchemes
 using DataFrames
 using EcologicalNetworksDynamics
-using Statistics
 
 # Define global community parameters.
 S = 20 # Species richness.
@@ -16,11 +17,9 @@ n_replicates = 100 # Number of food web replicates for each parameter combinatio
 C_values = [0.05, 0.1, 0.2] # Connectance values.
 tol_C = 0.01 # Tolerance on connectance when generating foodweb with the nichemodel.
 tmax = 2_000 # Simulation length.
-verbose = false # No information printed during simulation.
+verbose = false # Do not show '@info' messages during simulation run.
 
 """
-    standardize_K(foodweb, K, αij)
-
 Standardize total carrying capacity `K` for the number of producers in the `foodweb`
 and interspecific competition among producers `αij`.
 """
@@ -30,10 +29,10 @@ function standardize_K(foodweb, K, αij)
 end
 
 # Main simulation loop.
-# Each thread writes in its own DataFrame. We merge them at the end of the loop.
-df_to_merge = DataFrame[]
+# Each thread writes in its own DataFrame. Merge them at the end of the loop.
+dfs = [DataFrame() for _ in 1:length(C_values)] # Fill the vector with empty DataFrames.
 Threads.@threads for C in C_values # Parallelize on connctance values.
-    df_tmp = DataFrame(; C = Float64[], αij = Float64[], persistence = Float64[])
+    df_thread = DataFrame(; C = Float64[], αij = Float64[], persistence = Float64[])
     for i in 1:n_replicates
         foodweb = FoodWeb(nichemodel, S; Z, C, tol_C)
         for αij in αij_values
@@ -45,22 +44,22 @@ Threads.@threads for C in C_values # Parallelize on connctance values.
             # Measure species persistence i.e. the number of species that have
             # a biomass above 0 (`threshold`) at the last timestep (`last`).
             persistence = species_persistence(solution; threshold = 0, last = 1)
-            push!(df_tmp, [C, αij, persistence])
+            push!(df_thread, [C, αij, persistence])
             @info "C = $C, foodweb = $i, αij = $αij: done."
         end
     end
-    push!(df_to_merge, df_tmp)
+    push!(dfs, df_thread)
 end
-df = reduce(vcat, df_to_merge)
 @info "All simulations done."
+df = reduce(vcat, dfs)
 
-# We compute the mean persistence and the 95% confidence interval (`ic95`)
+# Compute the mean persistence and the 95% confidence interval (`ci95`)
 # for each (C, αij) combination.
 groups = groupby(df, [:C, :αij])
 df_processed = combine(
     groups,
     :persistence => mean,
-    :persistence => (x -> 1.96 * std(x) / sqrt(length(x))) => :ic95,
+    :persistence => (x -> 1.96 * std(x) / sqrt(length(x))) => :ci95,
 )
 
 # Plot mean species persistence with its confidence interval
@@ -73,17 +72,17 @@ ax = Axis(
     ylabel = "Species persistence",
 )
 curves = []
-colors = [get(ColorSchemes.viridis, val) for val in LinRange(0, 1, length(C_values))]
+colors = [get(viridis, val) for val in LinRange(0, 1, length(C_values))]
 for (C, color) in zip(C_values, colors)
     df_extract = df_processed[df_processed.C.==C, :]
     x = df_extract.αij
     y = df_extract.persistence_mean
-    ic = df_extract.ic95
+    ci = df_extract.ci95
     sl = scatterlines!(x, y; color = color, markercolor = color)
-    errorbars!(x, y, ic; color, whiskerwidth = 5)
+    errorbars!(x, y, ci; color, whiskerwidth = 5)
     push!(curves, sl)
 end
-Legend( # Put legend on top of the figure.
+Legend(
     fig[1, 1],
     curves,
     ["C = $C" for C in C_values];

use_cases/tritrophic-biomass-vs-temperature.jl

@@ -2,8 +2,8 @@ import AlgebraOfGraphics: set_aog_theme!
 import ColorSchemes: tol_light
 
 using CairoMakie
-using EcologicalNetworksDynamics
 using DataFrames
+using EcologicalNetworksDynamics
 
 # Create tri-tropic network.
 tritrophic = FoodWeb([0 0 0; 1 0 0; 0 1 0])
@@ -20,17 +20,18 @@ T_values = values(T0:1:T40)
 
 # DataFrame to store the simulations outputs.
 df = DataFrame(; T = Float64[], trophic_level = Int64[], Beq = Float64[])
+verbose = false # Do not show '@info' messages during simulation run.
 
 # Run simulations for each temperature across gradient.
 for T in T_values
     set_temperature!(params, T, ExponentialBA()) # Modify parameters with temperature.
     # Simulate biomass dynamics with modified parameters.
     B0 = params.environment.K[1] / 8 # Inital biomass.
     # `set_temperature!()` expresses rates in seconds,
-    # thus we also simulation time in seconds.
+    # then simulation time is expressed in seconds as well.
     # Simulation length and time steps are taken from Binzer et al., 2012.
     tmax = 315_360_000_000 # From Binzer et al., 2012.
-    callback = ExtinctionCallback(1e-6, false) # Remove TerminateSteadyState callback.
+    callback = ExtinctionCallback(1e-6, verbose) # Remove TerminateSteadyState callback.
     solution = simulate(params, [B0]; tmax, callback)
     Beq_vec = solution[end]
     for (Beq, tlvl) in zip(Beq_vec, trophic_lvl)
@@ -51,7 +52,7 @@ for (sub_df, color, markersize) in zip(groupby(df, :trophic_level), colors, mark
     sl = scatterlines!(T_celsius, sub_df.Beq; color, markercolor, markersize)
     push!(curves, sl)
 end
-Legend( # Put legend on top of the figure.
+Legend(
     fig[1, 1],
     curves,
     ["producer", "intermediate\nconsumer", "top predator"];