From de2bbcf91734f52091e88bf286b0cf8528c5b644 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Isma=C3=ABl=20Lajaaiti?= Date: Tue, 25 Apr 2023 18:02:09 +0200 Subject: [PATCH] Polishing Nutrient Intake. --- .vscode/settings.json | 2 - docs/src/example/paradox_enrichment.md | 4 +- docs/src/man/modelparameters.md | 76 ++- src/EcologicalNetworksDynamics.jl | 22 +- src/inputs/environment.jl | 34 +- src/inputs/producer_competition.jl | 107 ---- src/inputs/producer_growth.jl | 506 +++++++++--------- src/measures/functioning.jl | 25 +- src/measures/structure.jl | 131 +++-- src/model/dbdt.jl | 47 -- src/model/dudt.jl | 39 ++ src/model/generate_dbdt.jl | 49 +- src/model/model_parameters.jl | 45 +- src/model/nutrient_dynamics.jl | 18 - src/model/producer_growth.jl | 109 ++++ src/model/productivity.jl | 82 --- src/model/set_temperature.jl | 21 +- src/model/simulate.jl | 152 +++--- test/inputs/test-environment.jl | 9 +- test/inputs/test-producer_competition.jl | 52 +- test/measures/test-functioning.jl | 55 +- test/model/test-model_parameters.jl | 24 +- test/model/test-nutrient_intake.jl | 132 +++++ test/model/test-nutrientintake.jl | 114 ---- test/model/test-productivity.jl | 233 +++----- test/model/test-set_temperature.jl | 23 +- test/model/test-simulate.jl | 10 +- test/model/test-zombies.jl | 2 +- test/runtests.jl | 42 +- use_cases/diversity-vs-nti-intensity.jl | 8 +- use_cases/nutrients-competition-exclusion.jl | 65 +++ use_cases/paradox-enrichment.jl | 7 +- .../paradox_enrichment_nutrientintake.jl | 66 ++- .../persistence-vs-producer-competition.jl | 15 +- .../tritrophic-biomass-vs-temperature.jl | 7 +- use_cases/trophic-classes-vs-mortality.jl | 12 +- 36 files changed, 1160 insertions(+), 1185 deletions(-) delete mode 100644 .vscode/settings.json delete mode 100644 src/inputs/producer_competition.jl delete mode 100644 src/model/dbdt.jl create mode 100644 src/model/dudt.jl delete mode 100644 src/model/nutrient_dynamics.jl create mode 100644 src/model/producer_growth.jl delete mode 100644 src/model/productivity.jl create mode 100644 test/model/test-nutrient_intake.jl delete mode 100644 test/model/test-nutrientintake.jl create mode 100644 use_cases/nutrients-competition-exclusion.jl diff --git a/.vscode/settings.json b/.vscode/settings.json deleted file mode 100644 index 7a73a41b..00000000 --- a/.vscode/settings.json +++ /dev/null @@ -1,2 +0,0 @@ -{ -} \ No newline at end of file diff --git a/docs/src/example/paradox_enrichment.md b/docs/src/example/paradox_enrichment.md index a2b648e6..37f477b7 100644 --- a/docs/src/example/paradox_enrichment.md +++ b/docs/src/example/paradox_enrichment.md @@ -174,8 +174,8 @@ df = DataFrame(; # Run simulations: for each carrying capacity we compute the equilibrium biomass. for K in K_values - environment = Environment(foodweb; K) - params = ModelParameters(foodweb; functional_response, environment) + producer_growth = LogisticGrowth(foodweb; K) + params = ModelParameters(foodweb; functional_response, producer_growth) B0 = rand(2) # Inital biomass. solution = simulate(params, B0; tmax, verbose) extrema = biomass_extrema(solution, "10%") diff --git a/docs/src/man/modelparameters.md b/docs/src/man/modelparameters.md index 0bf64211..12aedf58 100644 --- a/docs/src/man/modelparameters.md +++ b/docs/src/man/modelparameters.md @@ -3,7 +3,7 @@ Once the [`FoodWeb`](@ref) is created, you still have to attribute values to the system parameters. To navigate easily through the parameters, -they are split in 4 different fields depending on their nature: +they are split into 4 different fields depending on their nature: - [`BioRates`](@ref) contains species biological rates, *e.g.* producer intrinsic growth rates @@ -11,7 +11,7 @@ they are split in 4 different fields depending on their nature: *e.g.* temperature - [`FunctionalResponse`](@ref) contains functional response parameters, *e.g.* the hill exponent - - [`producer_growth`](@ref) contains producers growth model parameters, + - [`producer_growth`](@ref) contains producer growth model parameters, *e.g.* the carrying capacity All the fields are stored within the same object [`ModelParameters`](@ref). @@ -123,12 +123,12 @@ can be given to [`BioRates`](@ref) as seen previously. biorates = BioRates(foodweb; x = rate) # custom metabolic demand ``` -## Producers growth +## Producer growth Two models are available for producers growth: -- a logistic growth model, used by default and a -- a nutrient intake model. + - a logistic growth model, used by default, and a + - a nutrient intake model. In both models, the producer biomasses $B_p$ grow following the general equation: @@ -136,21 +136,36 @@ $$ growth = r_i * B_i * G_i $$ -where $r_i$ is producer's $i$ intrinsic growth rate, $B_i$ its biomass and $G_i$ its net growth term. The latter can take different forms depending on the growth model in use. +Where $r_i$ is producer $i$'s intrinsic growth rate, $B_i$ is its biomass +and $G_i$ is its net growth term. +The latter can take different forms depending on the growth model in use. -The model is controlled by the field `producer_growth` in model parameters. +The form of producer growth term is controlled by +the field `producer_growth` in model parameters. ### The logistic model -In the logistic model, $G_i$ takes the general form +In the logistic model, $G_i$ takes the general form: $$ G_i = 1 - \frac{s}{K_i} $$ -where the numerator $s$ can take the simplest form $s = B_i$ (default) or include a matrix $\alpha$ describing the relative strength of the inter- vs intra-specific competition among producer. In this case, $s = \sum\alpha_{ij}B_j$ The matrix $\alpha$ should have as many rows and columns as there are producers in the food web. When the intra-specific competition is $1$ ($\alpha_{ii} = 1$), then species pairing with values below $1$ ($\alpha_{ij} < 1$) describe something akin to facilitation and values above $1$ ($\alpha_{ij} > 1$) describe a stronger inter-specific competition (that would lead to competitive exclusion in the absence of other processes). Note that when all non diagonal values are set to $0$, then $s = \sum\alpha_{ij}B_j = B_j$ and the two models are equivalent. +Where the numerator $s$ can take the simplest form $s = B_i$ (default), +or include a matrix $a$ describing the relative strength of the inter- vs +intra-specific competition among producers. +In this case, $s = \sum a_{ij}B_j$. -When setting up the simulations, the logistic model is controlled as followed: +The matrix $a$ should have as many rows +and columns as there are producers in the food web. +When the intra-specific competition is $1$ ($a_{ii} = 1$), +then species pairs with values below $1$ ($a_{ij} < 1$) describes something akin to +facilitation and values above $1$ ($a_{ij} > 1$) describe stronger inter-specific +competition (that would lead to competitive exclusion in the absence of other processes). +Note that when all non-diagonal values are set to $0$, +then $s = \sum\alpha_{ij}B_j = B_j$ and the two models are equivalent. + +When setting up the simulations, the logistic model is controlled as followed: ``` julia> foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]) @@ -161,10 +176,10 @@ FoodWeb of 3 species: method: unspecified species: [s1, s2, s3] -julia> logmod = LogisticGrowth(foodweb) #default behavior +julia> logmod = LogisticGrowth(foodweb) # Default behaviour. LogisticGrowth: Kᵢ - carrying capacity: [1.0, 1.0, nothing] - α - competition: (3, 3) sparse matrix + a - competition: (3, 3) sparse matrix julia> #Note: This is equivalent to ModelParameters(foodweb): julia> p = ModelParameters(foodweb, producer_growth = logmod) @@ -177,34 +192,44 @@ ModelParameters{BioenergeticResponse, LogisticGrowth}: temperature_response: NoTemperatureResponse ``` -which allows to control the value of the carrying capacities and the matrix of competition: +which allows for control over the value of the carrying capacities +and the matrix of competition: ``` julia> logmod = LogisticGrowth(foodweb, K = [10, 5, nothing], αij = 0.8) LogisticGrowth: Kᵢ - carrying capacity: [10.0, 5.0, nothing] - α - competition: (3, 3) sparse matrix + a - competition: (3, 3) sparse matrix ``` -### The nutrient intake model +### The nutrient intake model -The nutrient intake model is implemented following Brose et al., 2005. In this model, the net growth terms depends on nutrients $l$ concentration ($N_l$) and take the form: +The nutrient intake model is implemented following Brose et al., 2005. +In this model, the net growth terms depends on nutrient $l$'s concentration ($N_l$) +and take the form: $$ G_i(N) = MIN(\frac{N_l}{K_{li} + N_l}, ...) $$ -where $K_{li}$ describes the species-specific half saturation densities for each nutrient $l$. In the absence of other factors it describes the producers hierarchy of competition (lower values means higher intake efficiency). +where $K_{li}$ describes the species-specific half saturation densities +for each nutrient $l$. +In the absence of other factors it describes the producers hierarchy of competition +(lower values means higher intake efficiency). -The nutrients concentrations are described as: +The nutrients concentrations are described as: $$ dN_l/dt = D(S_l - N_l)-\sum^n_{i = 1}(C_{li}G_i(N)B_i) $$ -where $D$ is the system turnover (default is 0.25), $S_l$ is nutrient $l$'s supply concentration and $C_{li}$ is the concentration of the nutrients in the producers biomass (the higher it is the most needed the nutrient is for the species growth). +where $D$ is the system turnover (default is 0.25), $S_l$ is +nutrient $l$'s supply concentration and $C_{li}$ is the concentration of +the nutrients in the producers biomass +(the higher it is the more essential the nutrient is for the species growth). -As for the logistic model, the nutrient intake model's behavior is controlled through the use of the `producer_growth` field: +As with the logistic model, the nutrient intake model's behaviour is controlled +through the use of the `producer_growth` field: ``` julia> foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]) @@ -223,7 +248,8 @@ NutrientIntake - 2 nutrients: Kₗᵢ - half sat. densities: (S, n) matrix ``` -By default, the models has 2 nutrients, but has for every other parameters, this can be changed: +By default, the models has 2 nutrients, but as for every other parameters, +this can be changed: ``` julia> ni4 = NutrientIntake(foodweb, n = 4) @@ -234,7 +260,7 @@ NutrientIntake - 4 nutrients: Kₗᵢ - half sat. densities: (S, n) matrix ``` -It's then passed to `ModelParameter`: +This is then passed to `ModelParameter`: ``` julia> p = ModelParameters(foodweb, producer_growth = ni2) @@ -256,13 +282,13 @@ in which case the parameters take default values. ```@example econetd foodweb = FoodWeb([0 0 0; 0 0 0; 0 1 0]; Z = 10) -environment = Environment(foodweb) +environment = Environment() ``` By default the temperature is set to 293.15 K. ```@example econetd -environment.K +environment.T ``` However, these default values can be changed @@ -271,7 +297,7 @@ For the temperature you just have to provide the scalar corresponding to the wanted temperature. ```@example econetd -environment = Environment(foodweb; T = 273.15); +environment = Environment(; T = 273.15); environment.T ``` diff --git a/src/EcologicalNetworksDynamics.jl b/src/EcologicalNetworksDynamics.jl index 4eec884b..e09ac745 100644 --- a/src/EcologicalNetworksDynamics.jl +++ b/src/EcologicalNetworksDynamics.jl @@ -24,7 +24,6 @@ by writing `?` in a Julia REPL. """ module EcologicalNetworksDynamics -# Dependencies import DifferentialEquations.Rodas4 import DifferentialEquations.SSRootfind import DifferentialEquations.Tsit5 @@ -39,10 +38,8 @@ using Statistics using Decimals using SciMLBase -# Convenience aliases. const Solution = SciMLBase.AbstractODESolution -# Include scripts include(joinpath(".", "macros.jl")) include(joinpath(".", "inputs/foodwebs.jl")) include(joinpath(".", "inputs/nontrophic_interactions.jl")) @@ -50,16 +47,13 @@ include(joinpath(".", "inputs/functional_response.jl")) include(joinpath(".", "inputs/biological_rates.jl")) include(joinpath(".", "inputs/environment.jl")) include(joinpath(".", "inputs/temperature_dependent_rates.jl")) -#include(joinpath(".", "inputs/producer_competition.jl")) -include(joinpath(".", "inputs/nutrient_intake.jl")) include(joinpath(".", "inputs/producer_growth.jl")) include(joinpath(".", "model/model_parameters.jl")) +include(joinpath(".", "model/producer_growth.jl")) include(joinpath(".", "model/set_temperature.jl")) -#include(joinpath(".", "model/productivity.jl")) include(joinpath(".", "model/consumption.jl")) include(joinpath(".", "model/metabolic_loss.jl")) -include(joinpath(".", "model/dbdt.jl")) -include(joinpath(".", "model/nutrient_dynamics.jl")) +include(joinpath(".", "model/dudt.jl")) include(joinpath(".", "model/generate_dbdt.jl")) include(joinpath(".", "model/simulate.jl")) include(joinpath(".", "model/effect_nti.jl")) @@ -70,7 +64,6 @@ include(joinpath(".", "measures/utils.jl")) include(joinpath(".", "utils.jl")) include(joinpath(".", "formatting.jl")) -# Export public functions export @check_between export @check_greater_than export @check_in @@ -100,7 +93,6 @@ export DefaultGrowthParams export DefaultMaxConsumptionParams export DefaultMetabolismParams export DefaultMortalityParams -export DefaultNIntakeParams export draw_asymmetric_links export draw_symmetric_links export efficiency @@ -129,8 +121,12 @@ export get_parameters export handling_time export homogeneous_preference export interaction_names +export is_boostable export is_success export is_terminated +export ispredator +export isprey +export isproducer export Layer export LinearResponse export living_species @@ -149,6 +145,8 @@ export NIntakeParams export nontrophic_adjacency_matrix export NonTrophicIntensity export NoTemperatureResponse +export nutrient_indices +export nutrient_richness export NutrientIntake export population_stability export potential_competition_links @@ -167,15 +165,17 @@ export shannon_diversity export simpson export simulate export species_cv +export species_indices export species_persistence export species_richness export synchrony export TemperatureResponse export top_predators export total_biomass +export total_richness export trophic_classes export trophic_levels export trophic_structure export weighted_mean_trophic_level -end \ No newline at end of file +end diff --git a/src/inputs/environment.jl b/src/inputs/environment.jl index b82ac1ec..79d42c62 100644 --- a/src/inputs/environment.jl +++ b/src/inputs/environment.jl @@ -1,31 +1,25 @@ -#### Type definition #### mutable struct Environment - T::Union{Int64,Float64} + T::Float64 end -#### end #### -#### Type display #### """ One line Environment display. """ function Base.show(io::IO, environment::Environment) T = environment.T - print(io, "Environment(" * "T=$(T)K)") + print(io, "Environment(T=$(T)K)") end """ Multiline Environment display. """ function Base.show(io::IO, ::MIME"text/plain", environment::Environment) - - # Display output println(io, "Environment:") print(io, " T: $(environment.T) Kelvin") end -#### end #### """ - Environment(foodweb, T=293.15) + Environment(; T=293.15) Create environmental parameters of the system. @@ -38,23 +32,15 @@ The environmental parameters are: # Examples ```jldoctest -julia> foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]); # species 1 & 2 producers, 3 consumer - -julia> environment = Environment(foodweb) # default behaviour -Environment: - T: 293.15 Kelvin +julia> e = Environment() # Default behaviour. + e.T == 293.15 # Kelvin +true - -julia> Environment(foodweb; T = 300.15).T # change temperature -Environment: - T: 300.15 Kelvin +julia> e = Environment(; T = 300.0) + e.T == 300.0 +true ``` See also [`ModelParameters`](@ref). """ -function Environment( - net::EcologicalNetwork; - T::Real = 293.15, -) where {Tp<:Real} - Environment(T) -end +Environment(; T = 293.15) = Environment(T) diff --git a/src/inputs/producer_competition.jl b/src/inputs/producer_competition.jl deleted file mode 100644 index 4baef597..00000000 --- a/src/inputs/producer_competition.jl +++ /dev/null @@ -1,107 +0,0 @@ -#### Type definition #### -struct ProducerCompetition - α::SparseMatrixCSC{Float64,Int64} -end -#### end #### - -#### Type display #### -""" -One line Environment display. -""" -function Base.show(io::IO, alpha::ProducerCompetition) - print(io, "ProducerCompetition($(size(alpha.α)) matrix)") -end - -""" -Multiline Environment display. -""" -function Base.show(io::IO, ::MIME"text/plain", alpha::ProducerCompetition) - - # Display output - println(io, "Producer competition:") - println(io, " α: $(size(alpha.α)) matrix") - -end -#### end #### - -""" - ProducerCompetition(network, α = nothing, αii = 1.0, αij = 0.0) - -Create producer competition matrix of the system. - -The parameters are: - - - α the competition matrix of dimensions S*S, S being the species number of the - network. Set to nothing by default. - - αii the intracompetition term for all species, 1.0 by default. - - αij the interspecific competition term for all species, 0.0 by default. - -By default, the producers compete only with themselves (i.e. αii = 1.0, αij = 0.0). -In the resulting α matrix, the element α[i,j] represents the per capita -effect of species j on species i. If α matrix is specified, it overrides -the αii and αij parameters. Moreover, all the αij coefficients should be 0 for -non producers. - -# Examples - -```jldoctest -julia> A = [0 0 0; 0 0 0; 0 0 1] -3×3 Matrix{Int64}: - 0 0 0 - 0 0 0 - 0 0 1 - -julia> foodweb = FoodWeb(A; quiet = true) -FoodWeb of 3 species: - A: sparse matrix with 1 links - M: [1.0, 1.0, 1.0] - metabolic_class: 2 producers, 1 invertebrates, 0 vertebrates - method: unspecified - species: [s1, s2, s3] - -julia> c = ProducerCompetition(foodweb; αii = 0.5, αij = 1.0) -Producer competition: - α: (3, 3) matrix - -julia> my_α = [0.5 1.0 0; 1.0 0.5 0; 0 0 0] -3×3 Matrix{Float64}: - 0.5 1.0 0.0 - 1.0 0.5 0.0 - 0.0 0.0 0.0 - -julia> myc = ProducerCompetition(foodweb; α = my_α) -Producer competition: - α: (3, 3) matrix - -julia> c.α == myc.α -true -``` - -See also [`ModelParameters`](@ref). -""" -function ProducerCompetition(network::EcologicalNetwork; α = nothing, αii = 1.0, αij = 0.0) - # Matrix initialisation - S = richness(network) - non_producer = filter(!isproducer, network) - - if isnothing(α) - # Put the diagonal elements to αii - α = fill(αij, S, S) - for i in 1:S - α[i, i] = αii - end - - # Put coefficients of non-producers to 0 - α[non_producer, :] .= 0 - α[:, non_producer] .= 0 - - else - # α should be a square matrix - @assert size(α, 1) == size(α, 2) == S - # α should be 0 for non producers - @assert all(α[non_producer, :] .== 0) - @assert all(α[:, non_producer] .== 0) - end - - ProducerCompetition(SparseMatrixCSC(α)) -end diff --git a/src/inputs/producer_growth.jl b/src/inputs/producer_growth.jl index c0068794..3eb86319 100644 --- a/src/inputs/producer_growth.jl +++ b/src/inputs/producer_growth.jl @@ -1,308 +1,316 @@ -#= -Models for producer growth -=# - -#### Type definition #### abstract type ProducerGrowth end -struct LogisticGrowth <: ProducerGrowth - α::SparseMatrixCSC{Float64,Int64} # relative inter to intra-specific competition between producers - Kᵢ::Vector{Union{Nothing,Float64}} +""" +`LogisticGrowth` constructor. + +# Fields + + - `a`: producer competition matrix + - `K`: vector of producer carrying capacities + +See also [`NutrientIntake`](@ref). +""" +mutable struct LogisticGrowth <: ProducerGrowth + a::AbstractMatrix{Float64} + K::Vector{Union{Nothing,Float64}} end -struct NutrientIntake <: ProducerGrowth - n::Int64 - D::Real - Sₗ::Vector{Float64} - Cₗᵢ::SparseMatrixCSC{Float64} - Kₗᵢ::Matrix{Union{Nothing,Float64}} +mutable struct NutrientIntake <: ProducerGrowth + turnover::Vector{Float64} + supply::Vector{Float64} + concentration::Matrix{Float64} + half_saturation::Matrix{Union{Nothing,Float64}} end -#### end #### -#### Type display #### +""" + length(n::NutrientIntake) + +Number of nutrients in the `n`utrient intake model. +""" +Base.length(n::NutrientIntake) = length(n.turnover) + """ One line display FunctionalResponse """ -Base.show(io::IO, model::ProducerGrowth) = print(io, "$(typeof(model))") +Base.show(io::IO, model::ProducerGrowth) = show(io, typeof(model)) """ Multiline LogisticGrowth display. """ function Base.show(io::IO, ::MIME"text/plain", model::LogisticGrowth) - S = size(model.α, 1) + S = size(model.a, 1) println(io, "LogisticGrowth:") - println(io, " Kᵢ - carrying capacity: " * vector_to_string(model.Kᵢ)) - print(io, " α - competition: ($S, $S) sparse matrix") + println(io, " K - carrying capacity: " * vector_to_string(model.K)) + print(io, " a - competition: ($S, $S) matrix") end """ Multiline NutrientIntake display. """ function Base.show(io::IO, ::MIME"text/plain", model::NutrientIntake) - println(io, "NutrientIntake - $(model.n) nutrients:") - println(io, " D - turnover rate: $(model.D)") - println(io, " Sₗ - supply concentrations: " * vector_to_string(model.Sₗ)) - println(io, " Cₗᵢ - relative contents: (S, n) sparse matrix") - print(io, " Kₗᵢ - half sat. densities: (S, n) matrix") + n_nutrients = length(model) + n_producers = size(model.supply, 1) + println(io, "NutrientIntake - $n_nutrients nutrients:") + println(io, " `turnover` rate: $(model.turnover)") + println(io, " `supply` concentrations: " * vector_to_string(model.supply)) + println(io, " relative `concentration`: ($n_producers, $n_nutrients) matrix") + print(io, " `half_saturation` densities: ($n_producers, $n_nutrients) matrix") end -#### end #### - -#### Methods for building #### """ - LogisticGrowth(network, α = nothing, αii = 1.0, αij = 0.0, K = 1.0) - -Creates the parameters needed for the logistic producer growth model, these parameters are: - - α the competition matrix of dimensions S*S, S being the species number of the network. Set to nothing by default. - - αii the intracompetition term for all species, 1.0 by default. - - αij the interspecific (relative to intraspecific) competition term for all species, 0.0 by default. - - Kᵢ the vector of carrying capacities, set to 1 for all producer by default, meaning that each producer - gets its own carrying capacity and that the more producer in the system, the larger the overall carrying capacity. - -By default, with intraspecific competition set to unity and interspecific competition set to 0, the competition matrix -does not affect the system and the model behaves as if where not here. + LogisticGrowth( + network::EcologicalNetwork; + a = nothing, + K = 1, + quiet = false, + ) + +Create the parameters for the logistic producer growth model. +In the end, the `LogisticGrowth` struct created stores a vector of carrying capacities `K` +and a competition matrix between the producers `a`. + +The carrying capacities can be specified via the `K` arguments, +by default they are all set to 1 for the producers and to `nothing` for the consumers. + +By default the competition matrix `a` is assumed to have diagonal elements equal to 1 +and 0 for all off-diagonal elements. +This default can be changed via the `a` argument. +You can either directly pass the interaction matrix, +a single `Number` that will be interpreted as the value of all diagonal elements, +a tuple of `Number` the first will be interpreted as the diagonal value and +the second as the off-diagonal values, or a named tuple of `Number`. +For the latter, the following aliases can be used to refer to the diagonal elements +'diagonal', 'diag', 'd', and for the off-diagonal elements 'offdiagonal', 'offdiag', +'o', 'rest', 'nondiagonal', 'nond'. # Examples + ```jldoctest -julia> foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]); # species 1 & 2 producers, 3 consumer - -julia> def = LogisticGrowth(foodweb) # default behavior -LogisticGrowth: - K - carrying capacity: [1, 1, nothing] - α - competition: (3, 3) sparse matrix - -julia> def.α -3×3 SparseMatrixCSC{Float64, Int64} with 2 stored entries: - 1.0 ⋅ ⋅ - ⋅ 1.0 ⋅ - ⋅ ⋅ ⋅ - -julia> def.Kᵢ -3-element Vector{Union{Nothing, Real}}: - 1 - 1 - nothing - -julia> K_system = 10 #set system carrying capacity to 10 -10 - -julia> Kᵢ = K_system/length(producers(foodweb)) -5.0 - -julia> LogisticGrowth(foodweb; K = Kᵢ).Kᵢ # change the default value for producer carrying capacity -3-element Vector{Union{Nothing, Real}}: - 5.0 - 5.0 - nothing - -julia> LogisticGrowth(foodweb; K = [1, 2, nothing]).Kᵢ # can also provide a vector -3-element Vector{Union{Nothing, Real}}: - 1 - 2 - nothing - -julia> LogisticGrowth(foodweb; αii = 0.5, αij = 1.0).α -3×3 SparseMatrixCSC{Float64, Int64} with 4 stored entries: - 0.5 1.0 ⋅ - 1.0 0.5 ⋅ - ⋅ ⋅ ⋅ - -julia> #pass a matrix -julia> my_α = [0.5 1.0 0; 1.0 0.5 0; 0 0 0] -3×3 Matrix{Float64}: - 0.5 1.0 0.0 - 1.0 0.5 0.0 - 0.0 0.0 0.0 - -julia> myc = LogisticGrowth(foodweb; α = my_α).α -3×3 SparseMatrixCSC{Float64, Int64} with 4 stored entries: - 0.5 1.0 ⋅ - 1.0 0.5 ⋅ - ⋅ ⋅ ⋅ +julia> foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]) # 1 & 2 producers and 3 consumer. + g = LogisticGrowth(foodweb) + g.a == [ + 1 0 0 + 0 1 0 + 0 0 0 + ] +true + +julia> g.K == [1, 1, nothing] +true + +julia> K_cst = 5.0 + g = LogisticGrowth(foodweb; K = K_cst) # Change the default K. + g.K == [K_cst, K_cst, nothing] +true + +julia> K_vec = [1, 2, nothing] + g = LogisticGrowth(foodweb; K = K_vec) # Can also provide a vector. + g.K == K_vec +true + +julia> g = LogisticGrowth(foodweb; a = (diag = 0.5, offdiag = 1.0)) + g.a == [ + 0.5 1.0 0.0 + 1.0 0.5 0.0 + 0.0 0.0 0.0 + ] +true + +julia> my_competition_matrix = [ + 0.5 1.0 0.0 + 1.0 0.5 0.0 + 0.0 0.0 0.0 + ] + g = LogisticGrowth(foodweb; a = my_competition_matrix, quiet = true) + g.a == my_competition_matrix +true ``` -See also [`ModelParameters`](@ref). +See also [`ModelParameters`](@ref) and [`NutrientIntake`](@ref). """ -function LogisticGrowth(network::EcologicalNetwork; α = nothing, αii = 1.0, αij = 0.0, K = 1) +function LogisticGrowth(network::EcologicalNetwork; a = nothing, K = 1, quiet = false) S = richness(network) - - # carrying capacity isa(K, AbstractVector) || (K = [isproducer(i, network) ? K : nothing for i in 1:S]) @check_equal_richness length(K) S - - # competition - non_producer = filter(!isproducer, network) - if isnothing(α) - # Put the diagonal elements to αii - α = fill(αij, S, S) - for i in 1:S - α[i, i] = αii - end - # Put coefficients of non-producers to 0 - α[non_producer, :] .= 0 - α[:, non_producer] .= 0 + # Construct matrix from user input. + # Initialize these variables from the argument `a`. + diag, offdiag, raw = nothing, nothing, nothing + if isnothing(a) + diag, offdiag = 1.0, 0.0 # Default. + elseif a isa Number + diag, offdiag = a, 0.0 + elseif a isa AbstractMatrix # Raw matrix provided. + @assert size(a) == (S, S) + raw = a + elseif a isa Tuple && eltype(a) <: Real && length(a) in [1, 2] + n = length(a) + if n == 1 + diag, offdiag = a[1], 0.0 + elseif n == 2 + diag, offdiag = a + end + elseif a isa NamedTuple && eltype(a) <: Real && length(a) in [1, 2] + diag_aliases = [:diag, :diagonal, :d] # For convenience. + for key in diag_aliases + if haskey(a, key) + if !isnothing(diag) + @error "Ambiguous names in provided tuple: $(keys(a))\ + which one refers to the diagonal value?" + end + diag = a[key] + end + end + offdiag_aliases = [:offdiagonal, :offdiag, :o, :rest, :nondiagonal, :nd] + for key in offdiag_aliases + if haskey(a, key) + if !isnothing(offdiag) + @error "Ambiguous names in provided tuple: $(keys(a))\ + which one refers to the non-diagonal values?" + end + offdiag = a[key] + end + end + if isnothing(diag) || isnothing(offdiag) + throw( + ArgumentError( + "Wrong name in provided named tuple names $(keys(a)). \ + The valid names to specify the diagonal elements of the `a` matrix are \ + $diag_aliases and for the off-diagonal elements $offdiag_aliases.", + ), + ) + end + else + @error "Invalid matrix argument: $(a) of type $(typeof(a))." + end + # The potential information to build the competition matrix has been processed. + # Build the matrix or read it from the raw input. + if isnothing(raw) + prods = producers(network) + # If the offdiagonal elements are null (nothing or zero), then the matrix + # has mostly null entries, in which case a sparse matrix is used. + a = (isnothing(offdiag) || offdiag == 0) ? spzeros(S, S) : zeros(S, S) + for i in prods, j in prods + a[i, j] = i == j ? diag : offdiag + end else - # α should be a square matrix - @assert size(α, 1) == size(α, 2) == S - # α should be 0 for non producers - @assert all(α[non_producer, :] .== 0) - @assert all(α[:, non_producer] .== 0) + @assert size(raw) == (S, S) + a = sparse(raw) end - # build - LogisticGrowth(α, K) + # Competition between producers can only occur between a pair of producers. + non_producer = filter(!isproducer, network) + if !isempty(non_producer) + @assert all(a[non_producer, :] .== 0) + @assert all(a[:, non_producer] .== 0) + end + LogisticGrowth(a, K) end """ -NutrientIntake(network, n = 2, D = 0.25, Sₗ = repeat([10.0], n), Cₗ = [range(1, 0.5, length = n);], Kₗᵢ = 1.0) - -Creates the parameters needed for the nutrient intake producer growth model, these parameters are: - - n the number of nutrient in the model - - D the nutrients turnover rate - - Sₗ the supply concentration for each nutrient - - Cₗ the relative contents of the nutrients in producers - - Kₗᵢ the half-saturation densities for nutrient / producer pair + NutrientIntake( + network::EcologicalNetwork; + n_nutrients = 2, + turnover = 0.25, + supply = fill(10.0, n_nutrients), + concentration = collect(LinRange(1, 0.5, n_nutrients)), + half_saturation = 1.0, + ) + +Create parameters for the nutrient intake model of producer growths. +These parameters are: + + - the nutrient `turnover` rate (vector of length `n_nutrients`) + - the `supply` concentration for each nutrient (vector of length `n_nutrients`) + - the relative `concentration` of each nutrients in each producers + (matrix of size (n_producers, n_nutrients)) + - the `half_saturation` densities for each nutrient-producer pair + (matrix of size (n_producers, n_nutrients)) # Examples -```jldoctest -julia> foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]); # species 1 & 2 producers, 3 consumer - -julia> def = NutrientIntake(foodweb) # default behavior -NutrientIntake - 2 nutrients: - D - turnover rate: 0.25 - Sₗ - supply concentrations: [10.0, 10.0] - Cₗᵢ - relative contents: (2, 2) sparse matrix - Kₗᵢ - half sat. densities: (2, 2) sparse matrix -julia> def.n +```jldoctest +julia> foodweb = FoodWeb([4 => [1, 2, 3]]); # 4 feeds on 1, 2 and 3. + ni = NutrientIntake(foodweb) # knights who say.. + length(ni) # 2 nutrients by default. 2 -julia> def.D -0.25 +julia> ni.turnover == [0.25, 0.25] +true -julia> def.Sₗ -2-element Vector{Float64}: - 10.0 - 10.0 +julia> ni.supply == [10, 10] +true -julia> def.Cₗᵢ -2×2 SparseMatrixCSC{Float64, Int64} with 4 stored entries: - 1.0 0.5 - 1.0 0.5 +julia> ni.concentration == [1 0.5; 1 0.5; 1 0.5] +true -julia> def.Kₗᵢ -3×2 Matrix{Union{Nothing, Real}}: - 1.0 1.0 - 1.0 1.0 - nothing nothing +julia> ni.half_saturation == [1 1; 1 1; 1 1] +true -julia> # change in the number of producer affects other rates -julia> NutrientIntake(foodweb; n = 4).Kₗᵢ -3×4 Matrix{Union{Nothing, Real}}: - 1.0 1.0 1.0 1.0 - 1.0 1.0 1.0 1.0 - nothing nothing nothing nothing +julia> ni = NutrientIntake(foodweb; turnover = 1.0) + ni.turnover == [1.0, 1.0] +true +julia> ni = NutrientIntake(foodweb; n_nutrients = 5) + length(ni) +5 ``` -See also [`ModelParameters`](@ref). - +See also [`ModelParameters`](@ref) and [`ProducerGrowth`](@ref). """ -function NutrientIntake(network::EcologicalNetwork; - n = 2, - D = 0.25, - Sₗ = repeat([10.0], n), - Cₗᵢ = [range(1, 0.5, length = n);], - Kₗᵢ = 1.0) - - np = length(producers(network)) - S = richness(network) - - # Sanity check turnover (should be in ]0, 1]) - if (D <= 0) || (D > 1) - throw(ArgumentError("Turnover rate (D) should be in ]0, 1]")) - end - - # Check size of supply concentration and convert if needed - if length(Sₗ) == 1 - Sₗ = repeat([Sₗ], n) - elseif (length(Sₗ) != n) - throw(ArgumentError("Sₗ should have length n")) +function NutrientIntake( + network::EcologicalNetwork; + n_nutrients = 2, + turnover = 0.25, + supply = fill(10.0, n_nutrients), + concentration = collect(LinRange(1, 0.5, n_nutrients)), + half_saturation = 1.0, +) + 0 < turnover <= 1 || throw(ArgumentError("`turnover` rate should be in ]0, 1], \ + not $turnover.")) + + length(supply) == 1 && (supply = fill(supply[1], n_nutrients)) + length(supply) == n_nutrients || throw( + ArgumentError("`supply` should be of length 1 or `n_nutrients` = $n_nutrients."), + ) + + n_producers = length(producers(network)) + c = size(concentration) + if c == () || c == (1,) + concentration = fill(concentration[1], n_producers, n_nutrients) + elseif c == (n_nutrients,) + concentration = repeat(concentration, 1, n_producers)' + elseif c != (n_producers, n_nutrients) + throw(ArgumentError("`concentration` should be either: a number, \ + a vector of length 1 or number of producers = $n_producers, \ + or a matrix of size \ + (number of producer = $n_producers, \ + number of nutrients = $n_nutrients) \ + not $concentration (of type $(typeof(concentration))).")) end - # Convert C to array if needed - sc = size(Cₗᵢ) - if sc == (np, n) - Cₗᵢ = sparse(Cₗᵢ) - elseif sc == (n,) - C = repeat(Cₗᵢ, np) - Cₗᵢ = sparse(reshape(C, n, np) |> transpose) - else - throw(ArgumentError("Cₗᵢ should be a vector of length n or a matrix with dimensions (number of producer, n)")) - end - - #Check K and convert if needed - sk = size(Kₗᵢ) - if sk == (S, n) - K = Matrix{Union{Nothing,Float64}}(nothing, S, n) - K[producers(network), :] = Kₗᵢ[producers(network), :] - Kₗᵢ = K - elseif sk == (np, n) - K = Matrix{Union{Nothing,Float64}}(nothing, S, n) - K[producers(network), :] = Kₗᵢ - Kₗᵢ = K - elseif sk == (np,) || sk == () - K = Matrix{Union{Nothing,Float64}}(nothing, S, n) - K[producers(network), :] .= Kₗᵢ - Kₗᵢ = K - else - throw(ArgumentError("Kₗᵢ should be a Number, a vector of length n or a matrix with dimensions (number of producer, n)")) + length(turnover) == 1 && (turnover = fill(turnover[1], n_nutrients)) + length(turnover) != n_nutrients && + throw(ArgumentError("`turnover` should either be a number, a vector of length 1 \ + or a vector of length `n_nutrients` = $n_nutrients \ + not $turnover (of type $(typeof(turnover))).")) + + h = size(half_saturation) + if h == () || h == (1,) + half_saturation = fill(half_saturation[1], n_producers, n_nutrients) + elseif h == (n_producers,) + half_saturation = repeat(half_saturation, 1, n_nutrients) + elseif h != (n_producers, n_nutrients) + throw( + ArgumentError("`half_saturation` should be either: a number, \ + a vector of length 1 or number of producers = $n_producers, \ + or a matrix of size \ + (number of producer = $n_producers, \ + number of nutrients = $n_nutrients) \ + not $half_saturation (of type $(typeof(half_saturation))."), + ) end - - NutrientIntake(n, D, Sₗ, Cₗᵢ, Kₗᵢ) - -end - -#function NutrientIntake -#end -#### end #### - -# Producer Growth model functors -""" - LogisticGrowth(B, i, j) -""" -function (F::LogisticGrowth)(i, B, r, network::MultiplexNetwork, N = nothing) - r = effect_facilitation(r, i, B, network) - s = sum(F.α[i, :] .* B) - logisticgrowth(B[i], r[i], F.Kᵢ[i], s) -end -function (F::LogisticGrowth)(i, B, r, network::FoodWeb, N = nothing) - s = sum(F.α[i, :] .* B) - logisticgrowth(B[i], r[i], F.Kᵢ[i], s) -end -function logisticgrowth(B, r, K, s) - isnothing(K) && return 0 - r * B * (1 - s / K) -end - -""" - NutrientIntake(B, i, j) -""" -function (F::NutrientIntake)(i, B, r, network::MultiplexNetwork, N) - r = effect_facilitation(r, i, B, network) - nutrientintake(B[i], r[i], F.Kₗᵢ[i,:], N) -end -function (F::NutrientIntake)(i, B, r, network::FoodWeb, N) - nutrientintake(B[i], r[i], F.Kₗᵢ[i,:], N) -end -function nutrientintake(B, r, K, N) - all(isnothing.(K)) && return 0 - G_li = N ./ (K .+ N) - #if N = 0 -> G = NaN, but we want 0 (no growth) - G_li[isnan.(G_li)] .= 0.0 - minG = minimum(G_li) - r * B * minG + NutrientIntake(turnover, supply, concentration, half_saturation) end diff --git a/src/measures/functioning.jl b/src/measures/functioning.jl index ab193644..e15aa425 100644 --- a/src/measures/functioning.jl +++ b/src/measures/functioning.jl @@ -270,24 +270,23 @@ julia> foodweb = FoodWeb([0 1 1; 0 0 0; 0 0 0]); """ function producer_growth(solution; kwargs...) parameters = get_parameters(solution) - producer_idxs = producers(parameters.network) - producer_species = parameters.network.species[producer_idxs] + g = parameters.producer_growth + isa(g, LogisticGrowth) || throw(ArgumentError("Producer growth only makes sense for \ + logistic growth, and cannot be calculated with `$(typeof(g))` values.")) + prods = producers(parameters.network) # Producer indexes. + n_prods = length(prods) - Kp = parameters.environment.K[producer_idxs] - rp = parameters.biorates.r[producer_idxs] - αp = parameters.producer_competition.α[producer_idxs, producer_idxs] - - #Extract the producer_species biomass over the last timesteps - measure_on = extract_last_timesteps(solution; idxs = producer_idxs, kwargs...) + # Extract the producer biomass over the last timesteps. + measure_on = extract_last_timesteps(solution; kwargs...) + n_time_steps = size(measure_on, 2) - growth = zeros(length(producer_idxs), size(measure_on, 2)) - for (i, α) in enumerate(eachrow(αp)), (j, B) in enumerate(eachcol(measure_on)) - s = sum(α .* B) - growth[i, j] = logisticgrowth(B[i], rp[i], Kp[i], s) + growth = zeros(n_prods, n_time_steps) + for (i, p) in enumerate(prods), (j, B) in enumerate(eachcol(measure_on)) + growth[i, j] = g(p, B, parameters) end ( - species = producer_species, + species = parameters.network.species[prods], mean = mean.(eachrow(growth)), std = std.(eachrow(growth)), all = growth, diff --git a/src/measures/structure.jl b/src/measures/structure.jl index b65f3787..3e13b5ef 100644 --- a/src/measures/structure.jl +++ b/src/measures/structure.jl @@ -1,11 +1,89 @@ # Quantify food webs structural properties. """ + richness(p::ModelParameters) + Number of species in the network. + +```jldoctest +julia> foodweb = FoodWeb([1 => 2]) # 1 eats 2. + model = ModelParameters(foodweb) + richness(model) +2 +``` + +See also [`nutrient_richness`](@ref) and [`total_richness`](@ref). """ +richness(p::ModelParameters) = richness(p.network) richness(net::EcologicalNetwork) = richness(get_trophic_adjacency(net)) richness(A::AbstractMatrix) = size(A, 1) +""" +Species indices of the given `network`. +""" +species_indices(p::ModelParameters) = 1:richness(p) + +""" +Nutrient indices of the given `network`. +""" +function nutrient_indices(p::ModelParameters) + n = nutrient_richness(p) + S = richness(p) # Species richness. + S .+ (1:n) +end + +""" + nutrient_richness(p::ModelParameters) + +Number of nutrients in the model `p`. + +```jldoctest +julia> foodweb = FoodWeb([1 => 2]) # 1 eats 2. + producer_growth = LogisticGrowth(foodweb) # No nutrients. + model = ModelParameters(foodweb; producer_growth) + nutrient_richness(model) +0 + +julia> producer_growth = NutrientIntake(foodweb; n_nutrients = 3) # Include nutrients. + model_with_nutrients = ModelParameters(foodweb; producer_growth) + nutrient_richness(model_with_nutrients) +3 +``` + +See also [`richness`](@ref) and [`total_richness`](@ref). +""" +function nutrient_richness(p::ModelParameters) + pg = p.producer_growth + isa(pg, NutrientIntake) ? length(pg) : 0 +end + +""" + total_richness(p::ModelParameters) + +Total richness of the system defined as the sum of the number of species +and the number of nutrients. +If the model `p` does not include nutrient dynamics for producer growth, +then `total_richness` is equivalent the species [`richness`](@ref). + +```jldoctest +julia> foodweb = FoodWeb([1 => 2]) # 1 eats 2. + producer_growth = LogisticGrowth(foodweb) + model = ModelParameters(foodweb; producer_growth) + total_richness(model) == richness(model) == 2 +true + +julia> producer_growth = NutrientIntake(foodweb; n_nutrients = 3) # Include nutrients. + model_with_nutrients = ModelParameters(foodweb; producer_growth) + S = richness(model_with_nutrients) + N = nutrient_richness(model_with_nutrients) + total_richness(model_with_nutrients) == S + N == 5 # 2 + 3 +true +``` + +See also [`richness`](@ref) and [`nutrient_richness`](@ref). +""" +total_richness(p::ModelParameters) = richness(p) + nutrient_richness(p) + """ Connectance of network: number of links / (number of species)^2 """ @@ -13,7 +91,6 @@ connectance(A::AbstractMatrix) = sum(A) / richness(A)^2 connectance(foodweb::FoodWeb) = connectance(foodweb.A) connectance(U::UnipartiteNetwork) = connectance(U.edges) -#### Overloading Base methods #### """ Filter species of the network (`net`) for which `f(species_index, net) = true`. """ @@ -23,9 +100,7 @@ Base.filter(f, net::EcologicalNetwork) = filter(i -> f(i, net), 1:richness(net)) Transform species of the network (`net`) by applying `f` to each species. """ Base.map(f, net::EcologicalNetwork) = map(i -> f(i, net), 1:richness(net)) -#### end #### -#### Find producers #### """ Is species `i` of the network (`net`) a producer? """ @@ -34,17 +109,15 @@ isproducer(i, net::FoodWeb) = isproducer(i, net.A) isproducer(i, net::MultiplexNetwork) = isproducer(i, net.layers[:trophic].A) """ -Return indexes of the producers of the given `network`. +Return indices of the producers of the given `network`. """ producers(net) = filter(isproducer, net) producers(A::AbstractMatrix) = (1:richness(A))[all(A .== 0; dims = 2)|>vec] -#### end #### -#### Find predators #### """ predators_of(i, network) -Return indexes of the predators of species `i` for the given `network`. +Return indices of the predators of species `i` for the given `network`. # Examples @@ -78,7 +151,7 @@ ispredator(i, net::FoodWeb) = ispredator(i, net.A) ispredator(i, net::MultiplexNetwork) = ispredator(i, net.layers[:trophic].A) """ -Return indexes of the predators of the given `network`. +Return indices of the predators of the given `network`. """ predators(net::EcologicalNetwork) = filter(ispredator, net) #### end #### @@ -87,7 +160,7 @@ predators(net::EcologicalNetwork) = filter(ispredator, net) """ preys_of(i, network) -Return indexes of the preys of species `i` for the given `network`. +Return indices of the preys of species `i` for the given `network`. # Examples @@ -117,7 +190,7 @@ isprey(i, net::FoodWeb) = isprey(i, net.A) isprey(i, net::MultiplexNetwork) = isprey(i, net.layers[:trophic].A) """ -Return indexes of the preys of the network (`net`). +Return indices of the preys of the network (`net`). """ preys(net::EcologicalNetwork) = filter(isprey, net) preys(A::AbstractSparseMatrix) = [i for i in 1:size(A, 1) if !isempty(A[:, i].nzval)] @@ -128,9 +201,7 @@ Do species `i` and `j` share at least one prey? share_prey(i, j, A::AdjacencyMatrix) = !isempty(intersect(preys_of(i, A), preys_of(j, A))) share_prey(i, j, net::FoodWeb) = share_prey(i, j, net.A) share_prey(i, j, net::MultiplexNetwork) = share_prey(i, j, net.layers[:trophic].A) -#### end #### -#### Find verterbrates & invertebrates #### """ Is species `i` an ectotherm vertebrate? """ @@ -142,17 +213,15 @@ Is species `i` an invertebrate? isinvertebrate(i, net::EcologicalNetwork) = net.metabolic_class[i] == "invertebrate" """ -Return indexes of the vertebrates of the network (`net`). +Return indices of the vertebrates of the network (`net`). """ vertebrates(net::EcologicalNetwork) = filter(isvertebrate, net) """ -Return indexes of the invertebrates of the network (`net`). +Return indices of the invertebrates of the network (`net`). """ invertebrates(net::EcologicalNetwork) = filter(isinvertebrate, net) -#### end #### -#### Number of resources #### """ Number of resources species `i` is feeding on. """ @@ -166,7 +235,6 @@ Return a vector where element i is the number of resource(s) of species i. function number_of_resource(net::EcologicalNetwork) [number_of_resource(i, net) for i in 1:richness(net)] end -#### end #### function _ftl(A::AbstractMatrix) if isa(A, AbstractMatrix{Int64}) @@ -206,7 +274,7 @@ trophic_levels(net::EcologicalNetwork) = trophic_levels(get_trophic_adjacency(ne """ top_predators(net::EcologicalNetwork) -Top predator indexes (species eaten by nobody) of the given `net`work. +Top predator indices (species eaten by nobody) of the given `net`work. ```jldoctest julia> foodweb = FoodWeb([1 => 2, 2 => 3]) @@ -345,23 +413,25 @@ function remove_species(A::AbstractMatrix, species_to_remove) A_simplified end -function massratio(obj::Union{ModelParameters,FoodWeb}) - - if isa(obj, ModelParameters) - M = obj.foodweb.M - A = obj.foodweb.A - else - M = obj.M - A = obj.A - end +""" + mass_ratio(p::ModelParameters) - Z = mean((M./M')[findall(A)]) +Mean predator-prey body mass ratio given the model `p`arameters. +""" +mass_ratio(p::ModelParameters) = mass_ratio(p.network) - return Z +""" + mass_ratio(network::EcologicalNetwork) +Mean predator-prey body mass ratio given the `network`. +""" +function mass_ratio(network::EcologicalNetwork) + A = get_trophic_adjacency(network) + M = network.M + mean((M./M')[findall(A)]) end -#### Extend method from Graphs.jl to FoodWeb and MultiplexNetwork #### +# Extend method from Graphs.jl to FoodWeb and MultiplexNetwork. function Graphs.is_cyclic(net::EcologicalNetwork) is_cyclic(SimpleDiGraph(get_trophic_adjacency(net))) end @@ -369,4 +439,3 @@ end function Graphs.is_connected(net::EcologicalNetwork) is_connected(SimpleDiGraph(get_trophic_adjacency(net))) end -#### end #### diff --git a/src/model/dbdt.jl b/src/model/dbdt.jl deleted file mode 100644 index 83bed57a..00000000 --- a/src/model/dbdt.jl +++ /dev/null @@ -1,47 +0,0 @@ -#= -Core functions of the model -=# - -function dBdt!(dB, B, p, t) - - params, extinct_sp = p # unpack input - - S = richness(params.network) - - # Separate biomass and nutrients - N = B[S+1:end] - B = B[1:S] - - # Set up - Unpack parameters - response_matrix = params.functional_response(B, params.network) - r = params.biorates.r # vector of intrinsic growth rates - network = params.network - G = repeat([0.0], Int(S)) - - # Compute ODE terms for each species - for i in 1:S - growth = params.producer_growth(i, B, r, network, N) - G[i] = growth - eating, being_eaten = consumption(i, B, params, response_matrix) - metabolism_loss = metabolic_loss(i, B, params) - natural_death = natural_death_loss(i, B, params) - net_growth_rate = growth + eating - metabolism_loss - net_growth_rate = effect_competition(net_growth_rate, i, B, network) - dB[i] = net_growth_rate - being_eaten - natural_death - end - - for j in S+1:S+length(N) - l = j-S - dB[j] = nutrient_dynamics(params.producer_growth, l, B, N, G) - end - - # Avoid zombie species by forcing extinct biomasses to zero. - # https://github.com/BecksLab/EcologicalNetworksDynamics.jl/issues/65 - for sp in keys(extinct_sp) - B[sp] = 0.0 - end - - #recat B and N - B = vcat(B,N) -end - diff --git a/src/model/dudt.jl b/src/model/dudt.jl new file mode 100644 index 00000000..6bec7344 --- /dev/null +++ b/src/model/dudt.jl @@ -0,0 +1,39 @@ +""" + dudt!(du, u, p, _) + +Compute the species and nutrient (when relevant) abundance derivatives `du`, +given the abundances `u` and the model `p`. +The last silent argument is the time at which is evaluated the derivatives +and is a requirement of DifferentialEquations. +""" +function dudt!(du, u, p, _) + params, extinct_sp = p + S = richness(params) + B = u[species_indices(params)] + response_matrix = params.functional_response(B, params.network) + network = params.network + growth = fill(0.0, S) # Vector of producer growths. + + # Compute species biomass dynamics. + for i in species_indices(params) + growth[i] = params.producer_growth(i, u, params) + eating, being_eaten = consumption(i, B, params, response_matrix) + metabolism_loss = metabolic_loss(i, B, params) + natural_death = natural_death_loss(i, B, params) + net_growth_rate = growth[i] + eating - metabolism_loss + net_growth_rate = effect_competition(net_growth_rate, i, B, network) + du[i] = net_growth_rate - being_eaten - natural_death + end + + # Compute nutrient abundance dynamics. + for (i_nutrient, i_u) in enumerate(nutrient_indices(params)) + n = u[i_u] + du[i_u] = nutrient_dynamics(params, B, i_nutrient, n, growth) + end + + # Avoid zombie species by forcing extinct biomasses to zero. + # https://github.com/BecksLab/EcologicalNetworksDynamics.jl/issues/65 + for sp in keys(extinct_sp) + u[sp] = 0.0 + end +end diff --git a/src/model/generate_dbdt.jl b/src/model/generate_dbdt.jl index b210b9e0..b9152dbc 100644 --- a/src/model/generate_dbdt.jl +++ b/src/model/generate_dbdt.jl @@ -162,10 +162,9 @@ function generate_dbdt(parms::ModelParameters, type) style = Symbol(type) # TEMP: Summary of working and convincingly tested implementations. - resp = typeof(parms.functional_response) - net = typeof(parms.network) + (resp, _, pg) = boostable_criteria(parms) function to_test() - @warn "Automatic generated :$style specialized code for $resp ($net) \ + @warn "Automatic generated :$style specialized code for $resp ($net, $pg) \ has not been rigorously tested yet.\n\ If you are using it for non-trivial simulation, \ please make sure that the resulting trajectories \ @@ -174,25 +173,14 @@ function generate_dbdt(parms::ModelParameters, type) If they are, then consider adding that simulation to the packages tests set \ so this warning can be removed in future upgrades." end - unimplemented() = throw("Automatic generated :$style specialized code for $resp ($net) \ - is not implemented yet.") + unimplemented() = throw("Automatic generated :$style specialized code + for $resp ($net, $pg) is not implemented yet.") ok() = nothing - #! format: off - Dict( - (FoodWeb , :raw , LinearResponse) => to_test, - (FoodWeb , :raw , ClassicResponse) => to_test, - (FoodWeb , :raw , BioenergeticResponse) => to_test, - (FoodWeb , :compact, LinearResponse) => to_test, - (FoodWeb , :compact, ClassicResponse) => to_test, - (FoodWeb , :compact, BioenergeticResponse) => to_test, - (MultiplexNetwork, :raw , LinearResponse) => unimplemented, - (MultiplexNetwork, :raw , ClassicResponse) => unimplemented, - (MultiplexNetwork, :raw , BioenergeticResponse) => unimplemented, - (MultiplexNetwork, :compact, LinearResponse) => unimplemented, - (MultiplexNetwork, :compact, ClassicResponse) => unimplemented, - (MultiplexNetwork, :compact, BioenergeticResponse) => unimplemented, - )[(net, style, resp)]() - #! format: on + if !is_boostable(parms, type) + unimplemented() + else + to_test() # Switch to ok() once it has been tested on longer simulations. + end if style == :raw xp, data = generate_dbdt_raw(parms) @@ -206,5 +194,24 @@ function generate_dbdt(parms::ModelParameters, type) end +""" + is_boostable(parms::ModelParameters, type) + +Return true if boost has been implemented for this parametrization +with the given boosting style. +""" +function is_boostable(parms::ModelParameters, type) + (resp, _, pg) = boostable_criteria(parms) + resp != MultiplexNetwork && pg != NutrientIntake +end + +# Extract necessary information to decide whether boosting is supported or not. +function boostable_criteria(parms) + resp = typeof(parms.functional_response) + net = typeof(parms.network) + pg = typeof(parms.producer_growth) + (resp, net, pg) +end + include("./generate_dbdt_compact.jl") include("./generate_dbdt_raw.jl") diff --git a/src/model/model_parameters.jl b/src/model/model_parameters.jl index 17659581..a84b9dca 100644 --- a/src/model/model_parameters.jl +++ b/src/model/model_parameters.jl @@ -68,10 +68,9 @@ The parameters are compartmented in different groups: # Examples ```jldoctest -julia> foodweb = FoodWeb([0 1; 0 0]); # create a simple foodweb - -julia> p = ModelParameters(foodweb) -ModelParameters{BioenergeticResponse}: +julia> foodweb = FoodWeb([0 1; 0 0]) + p = ModelParameters(foodweb) +ModelParameters{BioenergeticResponse, LogisticGrowth}: network: FoodWeb(S=2, L=1) environment: Environment(T=293.15K) biorates: BioRates(d, r, x, y, e) @@ -79,7 +78,7 @@ ModelParameters{BioenergeticResponse}: producer_growth: LogisticGrowth temperature_response: NoTemperatureResponse -julia> p.network # check that stored foodweb is the same as the one we provided +julia> p.network # Check that stored foodweb is the same as the one we provided. FoodWeb of 2 species: A: sparse matrix with 1 links M: [1.0, 1.0] @@ -87,18 +86,16 @@ FoodWeb of 2 species: method: unspecified species: [s1, s2] -julia> p.functional_response # default is bioenergetic +julia> p.functional_response # Default is bioenergetic. BioenergeticResponse: B0: [0.5, 0.5] c: [0.0, 0.0] h: 2.0 ω: (2, 2) sparse matrix -julia> classic_response = ClassicResponse(foodweb); # choose classic functional response - -julia> p = ModelParameters(foodweb; functional_response = classic_response); - -julia> p.functional_response # check that the functional response is now "classic" +julia> classic_response = ClassicResponse(foodweb) + p = ModelParameters(foodweb; functional_response = classic_response); + p.functional_response # Check that the functional response is now "classic". ClassicResponse: c: [0.0, 0.0] h: 2.0 @@ -110,32 +107,30 @@ ClassicResponse: [`ModelParameters`](@ref) can also be generated from a [`MultiplexNetwork`](@ref). ```jldoctest -julia> foodweb = FoodWeb([0 1; 0 0]); # create a simple foodweb - -julia> net = MultiplexNetwork(foodweb); # convert to Multiplex Network - -julia> p = ModelParameters(net; functional_response = ClassicResponse(net)) -ModelParameters{ClassicResponse}: +julia> foodweb = FoodWeb([0 1; 0 0]) + net = MultiplexNetwork(foodweb) + p = ModelParameters(net; functional_response = ClassicResponse(net)) +ModelParameters{ClassicResponse, LogisticGrowth}: network: MultiplexNetwork(S=2, Lt=1, Lc=0, Lf=0, Li=0, Lr=0) - environment: Environment(K=[nothing, 1], T=293.15K) + environment: Environment(T=293.15K) biorates: BioRates(d, r, x, y, e) functional_response: ClassicResponse - producer_competition: ProducerCompetition((2, 2) matrix) + producer_growth: LogisticGrowth temperature_response: NoTemperatureResponse ``` """ function ModelParameters( network::EcologicalNetwork; biorates::BioRates = BioRates(network), - environment::Environment = Environment(network), + environment::Environment = Environment(), functional_response::FunctionalResponse = BioenergeticResponse(network), producer_growth::ProducerGrowth = LogisticGrowth(network), - temperature_response::TemperatureResponse = NoTemperatureResponse() + temperature_response::TemperatureResponse = NoTemperatureResponse(), ) - if isa(network, MultiplexNetwork) & !(isa(functional_response, ClassicResponse)) + if isa(network, MultiplexNetwork) && !(isa(functional_response, ClassicResponse)) type_response = typeof(functional_response) - @warn "Non-trophic interactions aren't implented for '$type_response'. - Use a functional response of type 'ClassicResponse' instead." + @warn "Non-trophic interactions for `$type_response` are not supported. \ + Use a classical functional response instead: `$ClassicResponse`." end ModelParameters( network, @@ -143,6 +138,6 @@ function ModelParameters( environment, functional_response, producer_growth, - temperature_response + temperature_response, ) end diff --git a/src/model/nutrient_dynamics.jl b/src/model/nutrient_dynamics.jl deleted file mode 100644 index 2a3fbdcb..00000000 --- a/src/model/nutrient_dynamics.jl +++ /dev/null @@ -1,18 +0,0 @@ -""" -**Nutrient uptake** -""" -function nutrient_dynamics(growthparams::NutrientIntake, l, B, N, G) - D = growthparams.D - S = growthparams.Sₗ[l] - cli = growthparams.Cₗᵢ[l] - return D * (S - N[l]) - sum(cli .* G .* B) -end - -""" -**Nutrient uptake** -Methode for logistic growth rate -""" -function nutrient_dynamics(growthparams::LogisticGrowth, l, B, N, G) - return 0.0 -end - diff --git a/src/model/producer_growth.jl b/src/model/producer_growth.jl new file mode 100644 index 00000000..1d56fc88 --- /dev/null +++ b/src/model/producer_growth.jl @@ -0,0 +1,109 @@ +# Producer growth functors. +function (g::LogisticGrowth)(i, u, params::ModelParameters) + isnothing(g.K[i]) && return 0.0 # Species i is not a producer. + B = u[species_indices(params)] + network = params.network + r = params.biorates.r + r_i = isa(network, MultiplexNetwork) ? effect_facilitation(r[i], i, B, network) : r[i] + s = sum(g.a[i, :] .* B) + r_i * B[i] * (1 - s / g.K[i]) +end + +function (g::NutrientIntake)(i, u, params::ModelParameters) + isproducer(i, params.network) || return 0.0 + network = params.network + prods = producers(network) + j = findfirst(prods .== i) + B = u[species_indices(params)] + N = u[nutrient_indices(params)] + r = params.biorates.r + r_i = isa(network, MultiplexNetwork) ? effect_facilitation(r[i], i, B, network) : r[i] + growth(N, k) = (N, k) == (0, 0) ? 0 : N / (N + k) + growth_vec = growth.(N, g.half_saturation[j, :]) + r_i * B[i] * minimum(growth_vec) +end + +# Code generation version (raw) (↑ ↑ ↑ DUPLICATED FROM ABOVE ↑ ↑ ↑). +# (update together as long as the two coexist) +function logisticgrowth(i, parms::ModelParameters) + B_i = :(B[$i]) + + # Non-producers are dismissed right away. + r_i = parms.biorates.r[i] + K_i = parms.producer_growth.K[i] + (r_i == 0 || isnothing(K_i)) && return 0 + + # Only consider nonzero competitors. + a_i = parms.producer_growth.a[i, :] + competitors = findall(a_i .!= 0) + a_i = a_i[competitors] + C = xp_sum([:c, :a], (competitors, a_i), :(a * B[c])) + + :($r_i * $B_i * (1 - $C / $K_i)) +end + +# Code generation version (compact): +# Explain how to efficiently construct all values of growth. +# This code assumes that dB[i] has already been *initialized*. +function growth(parms::ModelParameters, ::Symbol) + + # Pre-calculate skips over non-primary producers. + S = richness(parms.network) + data = Dict{Symbol,Any}( + :primary_producers => [ + (i, r, K) for (i, (r, K)) in + enumerate(zip(parms.biorates.r, parms.producer_growth.K)) if + (r != 0 && !isnothing(K)) + ], + ) + + # Flatten the e matrix with the same 'ij' indexing principle, + # but for producers competition links 'ic' instead of predation links. + # This should change in upcoming refactoring of compact code generation + # when MultiplexNetwork is supported, + # because a generic way accross all layers is needed for future compatibility. + a = parms.producer_growth.a + is, cs = findnz(SparseMatrixCSC(a)) # TODO: used to be sparse: fix perfs when dense? + data[:a] = [a[i, c] for (i, c) in zip(is, cs)] + data[:producers_competition_links] = (is, cs) + # Scratchspace to recalculate the sums on every timestep. + # /!\ Reset during iteration in 'consumption'. + # This will change in future evolution of :compact generated code + data[:s] = zeros(S) + + code = [ + :( + for (ic, (i, c)) in enumerate(zip(producers_competition_links...)) + s[i] += a[ic] * B[c] + end + ), + :( + for (i, r_i, K_i) in primary_producers # (skips over null terms) + dB[i] += r_i * B[i] * (1 - s[i] / K_i) + end + ), + ] + + code, data +end + +""" + nutrient_dynamics(model::ModelParameters, B, i_nutrients, n, G) + +Compute the dynamics of the nutrient `i_nutrient` given its abundance `n`, +the species biomass `B` and the vector of species growths `G` and the model `p`. + +The nutrient dynamics is applicable only if `p` is of type `NutrientIntake`. +""" +function nutrient_dynamics(model::ModelParameters, B, i_nutrient, n, G) + p = model.producer_growth + isp = producers(model.network) + if isa(p, LogisticGrowth) + throw(ArgumentError("Nutrient dynamics cannot be computed for producer growth \ + of type `$LogisticGrowth`.")) + end + d = p.turnover[i_nutrient] + s = p.supply[i_nutrient] + c = p.concentration[:, i_nutrient] + d * (s - n) - sum(c .* G[isp] .* B[isp]) +end diff --git a/src/model/productivity.jl b/src/model/productivity.jl deleted file mode 100644 index bfc3643d..00000000 --- a/src/model/productivity.jl +++ /dev/null @@ -1,82 +0,0 @@ -#= -Productivity -=# - -function logisticgrowth(i, B, r, K, s, network::MultiplexNetwork) - r = effect_facilitation(r, i, B, network) - logisticgrowth(B[i], r, K, s) -end - -logisticgrowth(i, B, r, K, s, _::FoodWeb) = logisticgrowth(B[i], r, K, s) -logisticgrowth(i, B, r, K, ::FoodWeb) = logisticgrowth(B[i], r, K, B[i]) - -logisticgrowth(i, B, r, K, network::MultiplexNetwork) = - logisticgrowth(i, B, r, K, B[i], network::MultiplexNetwork) - -function logisticgrowth(B, r, K, s = B) - !isnothing(K) || return 0 - r * B * (1 - s / K) -end -# Code generation version (raw) (↑ ↑ ↑ DUPLICATED FROM ABOVE ↑ ↑ ↑). -# (update together as long as the two coexist) -function logisticgrowth(i, parms::ModelParameters) - B_i = :(B[$i]) - - # Non-producers are dismissed right away. - r_i = parms.biorates.r[i] - K_i = parms.environment.K[i] - (r_i == 0 || isnothing(K_i)) && return 0 - - # Only consider nonzero competitors. - α_i = parms.producer_competition.α[i, :] - competitors = findall(α_i .!= 0) - α_i = α_i[competitors] - C = xp_sum([:c, :α], (competitors, α_i), :(α * B[c])) - - :($r_i * $B_i * (1 - $C / $K_i)) -end - -# Code generation version (compact): -# Explain how to efficiently construct all values of growth. -# This code assumes that dB[i] has already been *initialized*. -function growth(parms::ModelParameters, ::Symbol) - - # Pre-calculate skips over non-primary producers. - S = richness(parms.network) - data = Dict{Symbol,Any}( - :primary_producers => [ - (i, r, K) for - (i, (r, K)) in enumerate(zip(parms.biorates.r, parms.environment.K)) if - (r != 0 && !isnothing(K)) - ], - ) - - # Flatten the e matrix with the same 'ij' indexing principle, - # but for producers competition links 'ic' instead of predation links. - # This should change in upcoming refactoring of compact code generation - # when MultiplexNetwork is supported, - # because a generic way accross all layers is needed for future compatibility. - α = parms.producer_competition.α - is, cs = findnz(α) - data[:α] = [α[i, c] for (i, c) in zip(is, cs)] - data[:producers_competition_links] = (is, cs) - # Scratchspace to recalculate the sums on every timestep. - # /!\ Reset during iteration in 'consumption'. - # This will change in future evolution of :compact generated code - data[:s] = zeros(S) - - code = [ - :( - for (ic, (i, c)) in enumerate(zip(producers_competition_links...)) - s[i] += α[ic] * B[c] - end - ), - :( - for (i, r_i, K_i) in primary_producers # (skips over null terms) - dB[i] += r_i * B[i] * (1 - s[i] / K_i) - end - ), - ] - - code, data -end diff --git a/src/model/set_temperature.jl b/src/model/set_temperature.jl index 82b48380..2777ad37 100644 --- a/src/model/set_temperature.jl +++ b/src/model/set_temperature.jl @@ -1,27 +1,26 @@ #### Functors for temperature dependence methods #### -# No Temperature Response Functor -function (F::NoTemperatureResponse)(params::ModelParameters, T) - # record temperature in env, even though it has no effect - params.environment.T = T -end +# Record temperature in Environment even though it has no effect. +(F::NoTemperatureResponse)(params::ModelParameters, T) = params.environment.T = T # Exponential Boltzmann Arrhenius Functor. function (F::ExponentialBA)(params::ModelParameters, T;) + isa(params.producer_growth, NutrientIntake) && throw( + ArgumentError( + "Temperature dependence is not compatible with nutrient intake dynamics. \ + Either deactivate temperature dependence or \ + switch producer growth to `$LogisticGrowth`.", + ), + ) net = params.network - ## change params within BioRates params.biorates.r = Vector{Float64}(exp_ba_vector_rate(net, T, F.r)) params.biorates.x = Vector{Float64}(exp_ba_vector_rate(net, T, F.x)) - ## change params within FunctionalResponse params.functional_response.hₜ = exp_ba_matrix_rate(net, T, F.hₜ) params.functional_response.aᵣ = exp_ba_matrix_rate(net, T, F.aᵣ) - ## change params within Environment - params.environment.K = exp_ba_vector_rate(net, T, F.K) + params.producer_growth.K = exp_ba_vector_rate(net, T, F.K) params.environment.T = T end -### setting the temperature - # The entry point for the user. """ set_temperature!(p::ModelParameters, T, F!::TemperatureResponse) diff --git a/src/model/simulate.jl b/src/model/simulate.jl index 761f403b..e759c127 100644 --- a/src/model/simulate.jl +++ b/src/model/simulate.jl @@ -3,7 +3,7 @@ simulate( params::ModelParameters, B0::AbstractVector; - N0::AbstractVector = [1.0] + N0::AbstractVector = nothing, alg = nothing, t0::Number = 0, tmax::Number = 500, @@ -61,76 +61,45 @@ but want to find directly the biomass at steady state see [`find_steady_state`]( ```jldoctest julia> foodweb = FoodWeb([0 0; 1 0]); # create the foodweb - -julia> br = BioRates(foodweb; d = 0); # set natural death rate to 0 - -julia> params = ModelParameters(foodweb; biorates = br); - -julia> B0 = [0.5, 0.5]; # set initial biomass - -julia> solution = simulate(params, B0); # run simulation - -julia> is_terminated(solution) # => a steady state has been found -true - -julia> solution[begin] == B0 # initial biomass equals initial biomass -true - -julia> isapprox(solution[end], [0.188, 0.219]; atol = 1e-2) # steady state biomass -true + br = BioRates(foodweb; d = 0); # set natural death rate to 0 + params = ModelParameters(foodweb; biorates = br); + B0 = [0.5, 0.5]; # set initial biomass + solution = simulate(params, B0); # run simulation + @assert is_terminated(solution) # => a steady state has been found + @assert solution[begin] == B0 # initial biomass equals initial biomass + @assert isapprox(solution[end], [0.188, 0.219]; atol = 1e-2) # steady state biomass julia> using DifferentialEquations; - -julia> solution_custom_alg = simulate(params, B0; alg = BS5()); - -julia> isapprox(solution_custom_alg[end], [0.188, 0.219]; atol = 1e-2) -true + solution_custom_alg = simulate(params, B0; alg = BS5()); + @assert isapprox(solution_custom_alg[end], [0.188, 0.219]; atol = 1e-2) julia> import Logging; # TODO: remove when warnings removed from `generate_dbdt`. - -julia> # generate specialized code (same simulation) - -julia> xpr, data = Logging.with_logger(Logging.NullLogger()) do + # generate specialized code (same simulation) + xpr, data = Logging.with_logger(Logging.NullLogger()) do generate_dbdt(params, :raw) end; - -julia> solution = simulate(params, B0; diff_code_data = (eval(xpr), data)); - -julia> is_terminated(solution) # the same result is obtained, more efficiently. -true - -julia> solution[begin] == B0 -true - -julia> isapprox(solution[end], [0.188, 0.219]; atol = 1e-2) # steady state biomass -true + solution = simulate(params, B0; diff_code_data = (eval(xpr), data)); + @assert is_terminated(solution) # the same result is obtained, more efficiently. + @assert solution[begin] == B0 + @assert isapprox(solution[end], [0.188, 0.219]; atol = 1e-2) # steady state biomass julia> # Same with alternate style. - -julia> xpr, data = Logging.with_logger(Logging.NullLogger()) do + xpr, data = Logging.with_logger(Logging.NullLogger()) do generate_dbdt(params, :compact) end; + solution = simulate(params, B0; diff_code_data = (eval(xpr), data)); + @assert is_terminated(solution) + @assert solution[begin] == B0 + @assert isapprox(solution[end], [0.188, 0.219]; atol = 1e-2) # steady state biomass -julia> solution = simulate(params, B0; diff_code_data = (eval(xpr), data)); - -julia> is_terminated(solution) -true - -julia> solution[begin] == B0 -true - -julia> isapprox(solution[end], [0.188, 0.219]; atol = 1e-2) # steady state biomass -true ``` By default, the extinction callback throw a message when a species goes extinct. ```julia julia> foodweb = FoodWeb([0 0; 1 0]); - -julia> params = ModelParameters(foodweb; biorates = BioRates(foodweb; d = 0)); - -julia> simulate(params, [0.5, 1e-12]; verbose = true); # default: a message is thrown + params = ModelParameters(foodweb; biorates = BioRates(foodweb; d = 0)); + simulate(params, [0.5, 1e-12]; verbose = true); # default: a message is thrown [ Info: Species [2] is exinct. t=0.12316364776188903 julia> simulate(params, [0.5, 1e-12]; verbose = true); # no message thrown @@ -139,8 +108,8 @@ julia> simulate(params, [0.5, 1e-12]; verbose = true); # no message thrown """ function simulate( params::ModelParameters, - B0::AbstractVector; - N0::AbstractVector = [1.0], + B0; + N0 = nothing, alg = nothing, t0::Number = 0, tmax::Number = 500, @@ -148,22 +117,24 @@ function simulate( verbose = true, callback = CallbackSet( TerminateSteadyState(1e-6, 1e-4), - ExtinctionCallback(extinction_threshold, length(B0), verbose), + ExtinctionCallback(extinction_threshold, params, verbose), ), - diff_code_data = (dBdt!, params), + diff_code_data = (dudt!, params), kwargs..., ) - - # Check for consistency and format input arguments - S = richness(params.network) - all(B0 .>= 0) || throw( - ArgumentError( - "Inital biomasses provided in 'B0' should be all non-negative." * - "You gave B0 = $B0 which contains negative value(s).", - ), - ) - length(B0) == 1 && (B0 = repeat(B0, S)) + # Interpret parameters and check them for consistency. + S = richness(params) + all(B0 .>= 0) || + throw(ArgumentError("The argument received B0 = $B0 contains negative value(s). \ + Initial biomasses should all be non-negative.")) + length(B0) == 1 && (B0 = fill(B0[1], S)) @check_equal_richness length(B0) S + if !isnothing(N0) + all(N0 .>= 0) || throw( + ArgumentError("The argument received N0 = $N0 contains negative value(s). \ + Initial nutrient abundances should all be non-negative."), + ) + end @check_lower_than t0 tmax if isa(extinction_threshold, AbstractVector) length(extinction_threshold) == S || throw( @@ -173,7 +144,7 @@ function simulate( ) end - # Define ODE problem and solve + # Handle boosted simulations. timespan = (t0, tmax) code, data = diff_code_data # Work around julia's world count: @@ -190,22 +161,31 @@ function simulate( fun = (args...) -> Base.invokelatest(code, args...) extinct_sp = Dict(i => 0.0 for (i, b) in enumerate(B0) if b == 0.0) - # Set N0 to 0 if growth model is not nutrient intake. + # Set initial nutrient abundances `N0` and initial species biomass `B0`. if isa(params.producer_growth, NutrientIntake) - n = params.producer_growth.n - length(N0) == 1 && (N0 = repeat(N0, n)) - B0 = vcat(B0, N0) + isnothing(N0) && throw( + ArgumentError("If producer growth is of type `$NutrientIntake`, \ + use the `N0` argument to provide initial nutrient abundances."), + ) + n = nutrient_richness(params) + length(N0) == 1 && (N0 = fill(N0[1], n)) + @assert length(N0) == n || throw( + ArgumentError("The model contains $n nutrients, \ + but $(length(N0)) initial nutrient abundances were provided."), + ) + u0 = vcat(B0, N0) else - N0 = 0.0 + u0 = B0 end p = (params = data, extinct_sp = extinct_sp, original_params = params) - problem = ODEProblem(fun, B0, timespan, p) + timespan = (t0, tmax) + problem = ODEProblem(fun, u0, timespan, p) sol = solve( problem, alg; callback = callback, - isoutofdomain = (u, p, t) -> any(x -> x < 0, u[1:S]), + isoutofdomain = (u, p, t) -> any(x -> x < 0, u[species_indices(params)]), kwargs..., ) sol @@ -224,41 +204,39 @@ or an `AbstractVector` of length species richness (one threshold per species). If `verbose = true` a message is printed when a species goes extinct, otherwise no message are printed. """ -function ExtinctionCallback(extinction_threshold, S, verbose::Bool) - +function ExtinctionCallback(extinction_threshold, p::ModelParameters, verbose::Bool) # The callback is triggered whenever # a non-extinct species biomass goes below the threshold. - # Use either adequate code based on `extinction_threshold` type. # This avoids that the type condition be checked on every timestep. + sp = species_indices(p) # Vector of species indexes. species_under_threshold = if isa(extinction_threshold, Number) - (u, t, integrator) -> any(u[1:S][u[1:S].>0] .< extinction_threshold) + (u, _, _) -> any(u[sp][u[sp].>0] .< extinction_threshold) else - (u, t, integrator) -> any(u[1:S][u[1:S].>0] .< extinction_threshold[u[1:S].>0]) + (u, _, _) -> any(u[sp][u[sp].>0] .< extinction_threshold[u[sp].>0]) end get_extinct_species = if isa(extinction_threshold, Number) - (u, _) -> Set(findall(x -> x < extinction_threshold, u[1:S])) + (u, sp) -> Set(findall(x -> x < extinction_threshold, u[sp])) else - (u, S) -> Set((1:S)[u[1:S]. Set(sp[u[sp]. t for sp in new_extinct_sp) merge!(integrator.p.extinct_sp, new_extinct_sp_dict) # update extinct species list # Info message (printed only if verbose = true). if verbose && !isempty(new_extinct_sp) - S, S_ext = length(integrator.u[1:S]), length(all_extinct_sp) + S, S_ext = length(integrator.u[sp]), length(all_extinct_sp) @info "Species $([new_extinct_sp...]) went extinct at time t = $t. \n" * "$S_ext out of $S species are extinct." end diff --git a/test/inputs/test-environment.jl b/test/inputs/test-environment.jl index da37ba20..0bc56810 100644 --- a/test/inputs/test-environment.jl +++ b/test/inputs/test-environment.jl @@ -1,12 +1,7 @@ @testset "Environment" begin foodweb = FoodWeb([0 0 0; 1 0 0; 0 1 0]) - environment = Environment(foodweb) # default - @test environment.K == [1, nothing, nothing] + environment = Environment() # Default. @test environment.T == 293.15 - environment = Environment(foodweb; K = 10) # change K for producers (homogeneous) - @test environment.K == [10, nothing, nothing] - environment = Environment(foodweb; K = [10, 1, nothing]) # change K for producers (vec) - @test environment.K == [10, 1, nothing] - environment = Environment(foodweb; T = 300) # increase temperature + environment = Environment(; T = 300) # Custom temperature. @test environment.T == 300 end diff --git a/test/inputs/test-producer_competition.jl b/test/inputs/test-producer_competition.jl index e15abf56..a2f3a3d9 100644 --- a/test/inputs/test-producer_competition.jl +++ b/test/inputs/test-producer_competition.jl @@ -1,26 +1,42 @@ -@testset "ProducerCompetition" begin +@testset "Build LogisticGrowth with producer competition" begin foodweb = FoodWeb([0 0 0; 0 0 0; 0 1 0]; quiet = true) - compet = LogisticGrowth(foodweb) # default - @test compet.α == [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.0] - compet = LogisticGrowth(foodweb; αii = 1.0, αij = 1.0) # Put intercompetition - # Competition terms are 0 for all αij involving non-producers - @test compet.α == [1.0 1.0 0.0; 1.0 1.0 0.0; 0.0 0.0 0.0] - # Test that custom matrices work - myα = [1.0 1.0 0.0; 1.0 1.0 0.0; 0.0 0.0 0.0] - compet = LogisticGrowth(foodweb; α = myα) - @test compet.α == [1.0 1.0 0.0; 1.0 1.0 0.0; 0.0 0.0 0.0] - # Should fail if non zero α for non-producers - @test_throws AssertionError("all(α[non_producer, :] .== 0)") LogisticGrowth( + + g = LogisticGrowth(foodweb) # Default behaviour. + @test g.a == [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 0.0] + + # Define competition matrix with a (named) tuple. + for a in [(1, 2), (d = 1, nd = 2), (offdiag = 2, diag = 1), (diagonal = 1, rest = 2)] + g = LogisticGrowth(foodweb; a) # Set inter-specific competition. + # Competition terms are 0 for all αij involving non-producers + @test g.a == [1.0 2.0 0.0; 2.0 1.0 0.0; 0.0 0.0 0.0] + end + + # If invalid name(s) in named tuple an error is thrown. + @test_throws ArgumentError LogisticGrowth(foodweb; a = (d = 1, wrong_name = 2)) + @test_throws ArgumentError LogisticGrowth(foodweb; a = (oops = 1, ooops = 2)) + + # Custom competition matrix. + my_competition_matrix = [1.0 1.0 0.0; 1.0 1.0 0.0; 0.0 0.0 0.0] + g = LogisticGrowth(foodweb; a = my_competition_matrix, quiet = true) + @test g.a == [1.0 1.0 0.0; 1.0 1.0 0.0; 0.0 0.0 0.0] + + # Should fail if non zero a for non-producers. + message = "all(a[non_producer, :] .== 0)" + @test_throws AssertionError(message) LogisticGrowth( foodweb, - α = [1.0 1.0 0.0; 1.0 1.0 0.0; 1.0 0.0 0.0], + a = [1.0 1.0 0.0; 1.0 1.0 0.0; 1.0 0.0 0.0], + quiet = true, ) - # Should fail if the α matrix dim does not match the size of the foodweb - @test_throws AssertionError("size(α, 1) == size(α, 2) == S") LogisticGrowth( + # Should fail if the `a` dimension does not match the size of the foodweb. + message = "size(a) == (S, S)" + @test_throws AssertionError(message) LogisticGrowth( foodweb, - α = [1.0 0.0; 1.0 0.0; 0.0 0.0], + a = [1.0 0.0; 1.0 0.0; 0.0 0.0], + quiet = true, ) - @test_throws AssertionError("size(α, 1) == size(α, 2) == S") LogisticGrowth( + @test_throws AssertionError(message) LogisticGrowth( foodweb, - α = [1.0 0.0; 1.0 0.0], + a = [1.0 0.0; 1.0 0.0], + quiet = true, ) end diff --git a/test/measures/test-functioning.jl b/test/measures/test-functioning.jl index 38c8b0a2..66691478 100644 --- a/test/measures/test-functioning.jl +++ b/test/measures/test-functioning.jl @@ -118,45 +118,35 @@ end end @testset "Producer growth rate" begin - - # Set up foodweb = FoodWeb([0 0; 0 0]; quiet = true) params = ModelParameters(foodweb; biorates = BioRates(foodweb; d = 0)) - - sim = simulates(params, [0, 0.5]) - normal_growth = producer_growth(sim; last = 1) - # If biomass equal to 0, growth rate equal to 0 + u0 = [0, 0.5] # Initial biomass. + sol = simulates(params, u0) + normal_growth = producer_growth(sol; last = 1) + # If biomass equal to 0, growth rate equal to 0. @test normal_growth.all[normal_growth.species.=="s1"][1][1] ≈ 0.0 - first_growth = producer_growth(sim; last = length(sim.t), quiet = true) - # First growth rate should be equal to the logisticgrowth with initial - # biomass - params = get_parameters(sim) - s, r, K = params.network.species, params.biorates.r, params.environment.K - - @test first_growth.all[first_growth.species.=="s2", :][1] == - EcologicalNetworksDynamics.logisticgrowth.(0.5, r[s.=="s2"], K[s.=="s2"])[1] - @test first_growth.all[first_growth.species.=="s2", :][1] == 0.25 + # First growth rate should be equal to the LogisticGrowth with initial biomass. + first_growth = producer_growth(sol; last = length(sol.t), quiet = true) + model = get_parameters(sol) + g = model.producer_growth # Producer growth function. + @test isa(g, LogisticGrowth) + @test first_growth.all[1, :][1] == g(1, u0, model) == 0 + @test first_growth.all[2, :][1] == g(2, u0, model) == 0.25 - # Growth rate should converge to 0 as B converges to K - @test first_growth.all[first_growth.species.=="s2", :][length(sim.t)] ≈ 0 atol = 10^-3 + # Growth rate should converge to 0 as B converges to K. + @test first_growth.all[2, :][length(sol.t)] ≈ 0 atol = 10^-3 - # Growth rate computed only for producers + # Growth rate is computed only for producers. foodweb = FoodWeb([1 0; 0 0]; quiet = true) params = ModelParameters(foodweb) - sim = simulates(params, [0.5, 0.5]) - normal_growth = producer_growth(sim; last = 1) - + sol = simulates(params, [0.5, 0.5]) + normal_growth = producer_growth(sol; last = 1) @test length(normal_growth.species) == 1 - - # Test structure: @test length(normal_growth.mean) == length(normal_growth.std) == 1 - end @testset "Total biomass, species persistence, Hill numbers" begin - - # Set up foodweb = FoodWeb([0 0; 0 0]; quiet = true) params = ModelParameters(foodweb; biorates = BioRates(foodweb; d = 0)) @@ -165,22 +155,24 @@ end sim_zero = simulates(params, [0, 0]; verbose = false) m0, m1, m2 = sim_zero, sim_one_sp, sim_two_sp - # Total biomass should converge to K - tmp_K = get_parameters(m1).environment.K[2] + # Total biomass should converge to K. + tmp_K = get_parameters(m1).producer_growth.K[2] @test biomass(m1; last = 1).total ≈ tmp_K rtol = 0.001 @test biomass(m1; last = 1).total ≈ tmp_K rtol = 0.001 bm_two_sp = biomass(m2; last = 1) @test bm_two_sp.total == sum(bm_two_sp.species) @test bm_two_sp.total ≈ 2 atol = 0.001 - # Species sub selection works + # Select sub-collection of species. @test biomass(m1; last = 1, idxs = [1]).species[1] ≈ biomass(m1; last = 1, idxs = [1]).total @test biomass(m0; last = 2, quiet = true).total ≈ 0.0 - # Species richness is the binary equivalent of total_biomass + # Species richness is the binary equivalent of total_biomass. @test richness(m0; last = 2, quiet = true) ≈ 0.0 - @test species_persistence(m0; last = 2, quiet = true) ≈ 0.0 # Weird but it is a feature + + # Weird but it is a feature + @test species_persistence(m0; last = 2, quiet = true) ≈ 0.0 @test richness(m2[:, end]) == richness(m2; last = 1) == 2 @@ -263,5 +255,4 @@ end evennan([1, 1]; threshold = 1), ] @test all(evennan_check) - end diff --git a/test/model/test-model_parameters.jl b/test/model/test-model_parameters.jl index 368f377f..6342a911 100644 --- a/test/model/test-model_parameters.jl +++ b/test/model/test-model_parameters.jl @@ -2,43 +2,43 @@ A = [0 0 0; 1 0 0; 1 1 0] foodweb = FoodWeb(A) - # Default + # Default. p = ModelParameters(foodweb) @test p.biorates.x == [0, 0.314, 0.314] @test p.biorates.r == [1, 0, 0] - @test p.environment.K == [1, nothing, nothing] + @test p.producer_growth.K == [1, nothing, nothing] @test p.network.A == sparse(A) @test typeof(p.functional_response) == BioenergeticResponse - # Custom biorates + # Custom biorates. p = ModelParameters(foodweb; biorates = BioRates(foodweb; x = 1)) @test p.biorates.x == [1, 1, 1] # changed @test p.biorates.r == [1, 0, 0] # unchanged - @test p.environment.K == [1, nothing, nothing] # unchanged + @test p.producer_growth.K == [1, nothing, nothing] # unchanged @test p.network.A == sparse(A) # unchanged @test typeof(p.functional_response) == BioenergeticResponse # unchanged - # Classic Functional Response + # Classic Functional Response. p = ModelParameters(foodweb; functional_response = ClassicResponse(foodweb)) @test p.biorates.x == [0, 0.314, 0.314] # unchanged @test p.biorates.r == [1, 0, 0] # unchanged - @test p.environment.K == [1, nothing, nothing] # unchanged + @test p.producer_growth.K == [1, nothing, nothing] # unchanged @test p.network.A == sparse(A) # unchanged @test typeof(p.functional_response) == ClassicResponse # changed - # Linear Functional Response + # Linear Functional Response. p = ModelParameters(foodweb; functional_response = LinearResponse(foodweb)) @test typeof(p.functional_response) == LinearResponse - # Warning if not ClassicResponse and MultiplexNetwork + # Warning if not ClassicResponse and MultiplexNetwork. multiplex_network = MultiplexNetwork(foodweb) lresp = LinearResponse(multiplex_network) bresp = BioenergeticResponse(multiplex_network) cresp = ClassicResponse(multiplex_network) - linmsg = "Non-trophic interactions aren't implented for 'LinearResponse'. - Use a functional response of type 'ClassicResponse' instead." - biomsg = "Non-trophic interactions aren't implented for 'BioenergeticResponse'. - Use a functional response of type 'ClassicResponse' instead." + linmsg = "Non-trophic interactions for `LinearResponse` are not supported. \ + Use a classical functional response instead: `ClassicResponse`." + biomsg = "Non-trophic interactions for `BioenergeticResponse` are not supported. \ + Use a classical functional response instead: `ClassicResponse`." @test_logs (:warn, linmsg) ModelParameters( multiplex_network, functional_response = lresp, diff --git a/test/model/test-nutrient_intake.jl b/test/model/test-nutrient_intake.jl new file mode 100644 index 00000000..99d6642a --- /dev/null +++ b/test/model/test-nutrient_intake.jl @@ -0,0 +1,132 @@ +@testset "NutrientIntake: build struct." begin + + foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]) + ni = NutrientIntake(foodweb) + + # Check default values. + @test length(ni) == 2 + @test ni.turnover == [0.25, 0.25] + @test ni.supply == [10.0, 10.0] + @test ni.concentration == sparse([1 0.5; 1 0.5]) + @test ni.half_saturation == [1.0 1; 1 1] + + # Pass custom values. + ni = NutrientIntake( + foodweb; + n_nutrients = 3, + half_saturation = 10, + turnover = 0.8, + supply = 4, + concentration = [0.4, 0.1, 0.5], + ) + @test length(ni) == 3 + @test ni.turnover == fill(0.8, 3) + @test ni.supply == repeat([4.0], length(ni)) + @test ni.concentration == sparse([0.4 0.1 0.5; 0.4 0.1 0.5]) + @test ni.half_saturation == [10.0 10.0 10.0; 10.0 10.0 10.0] + + ni_1 = NutrientIntake( + foodweb; + n_nutrients = 4, + half_saturation = [1, 1], + supply = [1, 1, 1, 1], + concentration = [0.2, 0.2, 0.2, 0.2], + ) + ni_2 = NutrientIntake( + foodweb; + n_nutrients = 4, + half_saturation = ones(2, 4), + supply = 1, + concentration = 0.2, + ) + @test length(ni_1) == length(ni_2) == 4 + @test ni_1.turnover == ni_2.turnover == fill(0.25, 4) + + # `half_saturation` is a matrix of size (n_producers, n_nutrients) + @test ni_1.half_saturation == ni_2.half_saturation == [ # size = (n_prod, n_nutrients). + 1.0 1.0 1.0 1.0 + 1.0 1.0 1.0 1.0 + ] + # `supply` is a vector of length n_nutrients. + @test ni_1.supply == ni_2.supply == [1.0, 1.0, 1.0, 1.0] + # `concentration` is a matrix of size (n_producers, n_nutrients) + @test ni_1.concentration == ni_2.concentration == fill(0.2, 2, 4) + + # Error if arguments have the wrong dimensions. + @test_throws ArgumentError NutrientIntake(foodweb, half_saturation = [1, 1, 1]) + @test_throws ArgumentError NutrientIntake(foodweb, half_saturation = [1 1 1; 1 1 1]) + @test_throws ArgumentError NutrientIntake(foodweb, supply = [1 1 1]) + @test_throws ArgumentError NutrientIntake(foodweb, concentration = [1 1 1]) + @test_throws ArgumentError NutrientIntake( + foodweb, + n_nutrients = 3, + concentration = [1 1; 1 1; 1 1], + ) +end + +@testset "NutrientIntake: functor." begin + + foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]) + ni = NutrientIntake(foodweb) + model = ModelParameters(foodweb; producer_growth = ni) + n = total_richness(model) + u = fill(0, n) # Species biomass and nutrient abundances. + # If all abundances are set to 0, then the growth term is null. + @test ni(1, u, model) == ni(2, u, model) == 0 + + # If nutrient abundance and species biomass are strictly positive, + # then the growth is also positive and increase the producer biomass. + N = ones(2) # Nutrient abundances. + B = [1, 2, rand()] # Species abundances. + u = vcat(B, N) + @test 0 < ni(1, u, model) < ni(2, u, model) + @test ni(3, u, model) == 0 # Growth term is null non-producers. + + # Test a specific case. + N = rand(2) + B = rand(3) + u = vcat(B, N) + half_saturation = rand(2) + ni = NutrientIntake(foodweb; half_saturation) + model = ModelParameters(foodweb; producer_growth = ni) + expectations = [B[i] * minimum(N ./ (ni.half_saturation[i, :] .+ N)) for i in 1:2] + push!(expectations, 0) # Zero growth is expected for species 3 (consumer). + for i in 1:3 + @test ni(i, u, model) == expectations[i] + end +end + +@testset "NutrientIntake: dynamics." begin + + # Consumer goes extinct if there is no supply. + # Producer do not as long as there are no metabolic losses or mortality outside + # of consumption. + foodweb = FoodWeb([0 0 0; 1 0 0; 0 1 0]) + ni = NutrientIntake(foodweb; supply = 0.0) + S = richness(foodweb) + model = ModelParameters(foodweb; producer_growth = ni) + B0 = ones(S) + N0 = ones(nutrient_richness(model)) + sim = simulates(model, B0; N0) + @test sim[end][2:3] == [0.0, 0.0] + + # In the absence of consumers, all biomasses are equal (for equal growth parameters). + foodweb = FoodWeb([0 0; 0 0]) + producer_growth = NutrientIntake(foodweb) + model = ModelParameters(foodweb; producer_growth) + B0 = ones(2) + N0 = ones(2) + sol = simulates(model, B0; N0) + traj = reduce(hcat, sol.u) # row = species & nutrients, col = time steps. + sp = species_indices(model) # Species indexes. + @test all(traj[sp[1], :] .== traj[sp[2], :]) + + # `half_saturation` sets the hierarchy of competition between producers. + # The smaller the half saturation, the larger the producer grows. + producer_growth = NutrientIntake(foodweb; half_saturation = [0.1 0.1; 0.5 0.5]) + model = ModelParameters(foodweb; producer_growth) + B0 = ones(2) + N0 = ones(2) + sol = simulates(model, B0; N0) + @test sol[end][1] > sol[end][2] +end diff --git a/test/model/test-nutrientintake.jl b/test/model/test-nutrientintake.jl deleted file mode 100644 index 7153fbed..00000000 --- a/test/model/test-nutrientintake.jl +++ /dev/null @@ -1,114 +0,0 @@ -@testset "Logistic growth" begin - - #= - OBJECT CONSTRUCTION - =# - - # default - foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]) - ni = NutrientIntake(foodweb) - - # check default values - @test ni.n == 2 #defult is 2 nutrients - @test ni.D == 0.25 - @test ni.Sₗ == [10.0, 10.0] - @test ni.Cₗᵢ == sparse([1 0.5 ; 1 0.5]) - @test ni.Kₗᵢ == [1.0 1 ; 1 1 ; nothing nothing] - - # passing custom values - n changes the dimensions of C, S, and K - ni = NutrientIntake(foodweb, n = 3, Kₗᵢ = 10, D = 0.8, Sₗ = 4, Cₗᵢ = [0.4, 0.1, 0.5]) - @test ni.n == 3 - @test ni.D == 0.8 - @test ni.Sₗ == repeat([4.0], ni.n) - @test ni.Cₗᵢ == sparse([0.4 0.1 0.5 ; 0.4 0.1 0.5]) - @test ni.Kₗᵢ == [10.0 10.0 10.0 ; 10.0 10.0 10.0 ; nothing nothing nothing] - - # passing custom values - ni1 = NutrientIntake(foodweb, n = 4, Kₗᵢ = [1, 1], Sₗ = [1, 1, 1, 1], Cₗᵢ = [0.2, 0.2, 0.2, 0.2]) - ni2 = NutrientIntake(foodweb, n = 4, Kₗᵢ = ones(2,4), Sₗ = ones(4), Cₗᵢ = fill(0.2, 2, 4)) - ni3 = NutrientIntake(foodweb, n = 4, Kₗᵢ = ones(3,4)) - @test ni1.n == ni2.n == ni3.n == 4 - @test ni1.D == ni2.D == ni3.D == 0.25 - # - K has dimenion #species, #nutrients (K == nothing for non prod.) - @test ni1.Kₗᵢ == ni2.Kₗᵢ == ni3.Kₗᵢ == - [1.0 1.0 1.0 1.0 ; - 1.0 1.0 1.0 1.0 ; - nothing nothing nothing nothing] - # - Sₗ has dimension #producer, #nutrient - @test ni1.Sₗ == ni2.Sₗ == - [1.0, 1.0, 1.0, 1.0] - # - Cₗᵢ has dimension #producer, #nutrient - @test ni1.Cₗᵢ == ni2.Cₗᵢ == - fill(0.2, 2, 4) - - # Error for wrong dimensions - @test_throws ArgumentError NutrientIntake(foodweb, Kₗᵢ = [1,1,1]) - @test_throws ArgumentError NutrientIntake(foodweb, Kₗᵢ = [1 1 1 ; 1 1 1]) - @test_throws ArgumentError NutrientIntake(foodweb, Sₗ = [1 1 1]) - @test_throws ArgumentError NutrientIntake(foodweb, Cₗᵢ = [1 1 1]) - @test_throws ArgumentError NutrientIntake(foodweb, n = 3, Cₗᵢ = [1 1 ; 1 1 ; 1 1]) - - #= - GROWTH FUNCTION - =# - - #If N = 0, growth is null - S = richness(foodweb) - B = ones(S) - r = BioRates(foodweb).r - N = zeros(2) - ni = NutrientIntake(foodweb) - - @test ni(1, B, r, foodweb, N) == 0.0 - @test ni(2, B, r, foodweb, N) == 0.0 - - # Sanity check: - #If N != 0, growth is positive for default values and dependent on B - N = ones(2) - B = [3,2,1] - @test ni(1, B, r, foodweb, N) > ni(2, B, r, foodweb, N) > 0.0 - # Still 0 for non producers though - @test ni(3, B, r, foodweb, N) == 0 - - # Test a specific case - N = [2.5, 3.2] - B = [0.3, 0.5, 0.2] - ni = NutrientIntake(foodweb, Kₗᵢ = [1.2, 3.6]) - exp1 = r[1] * B[1] * minimum(N ./ (ni.Kₗᵢ[1,:] .+ N)) - exp2 = r[2] * B[2] * minimum(N ./ (ni.Kₗᵢ[2,:] .+ N)) - exp3 = 0 - @test ni(1, B, r, foodweb, N) == exp1 - @test ni(2, B, r, foodweb, N) == exp2 - @test ni(3, B, r, foodweb, N) == exp3 - - #= - DYNAMICS - =# - - # Sanity check: - # consumers goes extinct if supply is null - # (producer don't because there are no metabolic losses or mortality outside of consumption) - foodweb = FoodWeb([0 0 0; 1 0 0; 0 1 0]) - ni = NutrientIntake(foodweb, Sₗ = 0.0) - S = richness(foodweb) - p = ModelParameters(foodweb, producer_growth = ni) - B0 = ones(S) - sim = simulate(p, B0) - - @test sim[end][2:3] == [0.0, 0.0] - - #in the absence of consumers, all biomasses are equal (for equal growth parameters) - foodweb = FoodWeb([0 0 ; 0 0]) - ni = NutrientIntake(foodweb) - p = ModelParameters(foodweb, producer_growth = ni) - sim = simulate(p, ones(2)) - - @test all(getindex.(sim.u, 1) .== getindex.(sim.u, 2)) - - #Kₗᵢ sets the hierarchy of competition between producers - ni = NutrientIntake(foodweb, Kₗᵢ = [0.1 0.1 ; 0.5 0.5]) #smaller value = better - p = ModelParameters(foodweb, producer_growth = ni) - sim = simulate(p, ones(2)) - - @test sim[end][1] > sim[end][2] -end \ No newline at end of file diff --git a/test/model/test-productivity.jl b/test/model/test-productivity.jl index c6dea4a8..338d13b1 100644 --- a/test/model/test-productivity.jl +++ b/test/model/test-productivity.jl @@ -1,181 +1,70 @@ -@testset "Logistic growth" begin - # Intern method - B, r, K, s = 1, 1, nothing, 0 - @test EcologicalNetworksDynamics.logisticgrowth(B, r, K, s) == 0 # K is nothing, growth is null - B, r, K, s = 1, 0, 1, 0 - @test EcologicalNetworksDynamics.logisticgrowth(B, r, K, s) == 0 # r is null, growth is null - B, r, K, s = 0, 1, 1, 0 - @test EcologicalNetworksDynamics.logisticgrowth(B, r, K, s) == 0 # B is null, growth is null - B, r, K, s = 2, 1, 2, 2 - @test EcologicalNetworksDynamics.logisticgrowth(B, r, K, s) == 0 # B = K, growth is null - @test EcologicalNetworksDynamics.logisticgrowth.(0.5, 1, 1, 0.5) == 0.5 * 1 * (1 - 0.5 / 1) +@testset "LogisticGrowth functor" begin - # Extern method without facilitation and with intracompetition only - foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]) # 1 & 2 producers - S = EcologicalNetworksDynamics.richness(foodweb) - p = ModelParameters( - foodweb; - producer_growth = LogisticGrowth(foodweb; αii = 1.0, αij = 0.0), - ) - K, r, α = p.environment.K, p.biorates.r, p.producer_growth.α + foodweb_to_test = [ + FoodWeb([0 0; 0 0]; quiet = true), # 2 producers. + FoodWeb([0 0 0; 0 0 0; 1 1 0]; quiet = true), # 2 producers and 1 consumer. + ] - mat_expected_growth = [0 0 0; 0.25 0.25 0] - for (i, B) in enumerate(fill.([1, 0.5], S)) # ~ [[1, 1, 1], [0.5, 0.5, 0.5]] - for (sp, expected_growth) in enumerate(mat_expected_growth[i, :]) - println((i, sp)) - prodgrowth = LogisticGrowth(foodweb) - prodgrowth_α = LogisticGrowth(foodweb, α = α) - @test prodgrowth(sp, B, r, foodweb, nothing) == - prodgrowth_α(sp, B, r, foodweb, nothing) == - expected_growth - end - end + for (i, network) in enumerate(foodweb_to_test), f0 in [0, rand()] + # Default behaviour. + S = richness(network) + A = zeros(Integer, S, S) + A[1, 2] = A[2, 1] = 1 + f0 > 0 && (network = MultiplexNetwork(network; facilitation = (A = A, I = f0))) + functional_response = ClassicResponse(network) + model = ModelParameters(network; functional_response) + g = model.producer_growth + @test isa(g, LogisticGrowth) + u = rand(richness(model)) + @test g(1, u, model) == (1 + f0 * u[2]) * u[1] * (1 - u[1] / 1) + @test g(2, u, model) == (1 + f0 * u[1]) * u[2] * (1 - u[2] / 1) + i == 2 && @test g(3, u, model) == 0.0 - p = ModelParameters(foodweb; biorates = BioRates(foodweb; r = 2)) - K, r, α = p.environment.K, p.biorates.r, p.producer_growth.α - B = [0.5, 0.5, 0.5] - expected_growth = [0.5, 0.5, 0] - for i in 1:3 - prodgrowth = LogisticGrowth(foodweb) - prodgrowth_α = LogisticGrowth(foodweb, α = α) - @test prodgrowth(i, B, r, foodweb, nothing) == - prodgrowth_α(i, B, r, foodweb, nothing) == - expected_growth[i] - end + # Change carrying capacity and intrinsic growth rate. + K = [isproducer(i, network) ? 1 + rand() : nothing for i in species_indices(model)] + r = [isproducer(i, network) ? rand() : 0 for i in species_indices(model)] + g = LogisticGrowth(network; K) + biorates = BioRates(network; r) + model = ModelParameters(network; producer_growth = g, biorates, functional_response) + @test g(1, u, model) == r[1] * (1 + f0 * u[2]) * u[1] * (1 - u[1] / K[1]) + @test g(2, u, model) == r[2] * (1 + f0 * u[1]) * u[2] * (1 - u[2] / K[2]) + i == 2 && @test g(3, u, model) == 0.0 - # Extern method with intra and intercompetition - p = ModelParameters( - foodweb; - producer_growth = LogisticGrowth(foodweb; αii = 1.0, αij = 1.0), - ) - K, r, α = p.environment.K, p.biorates.r, p.producer_growth.α - B = [0.5, 0.5, 0.5] - for i in 1:3 - # It is like each producer had a density of one (look the "B * 2" in the - # left call of the function) - prodgrowth = LogisticGrowth(foodweb) - prodgrowth_α = LogisticGrowth(foodweb, α = α) - @test prodgrowth(i, B*2, r, foodweb, nothing) == - prodgrowth_α(i, B, r, foodweb, nothing) == - 0.0 + # With producer competition. + # Change intra-specific competition only. + a_ii = rand() + g = LogisticGrowth(network; a = a_ii) + model = ModelParameters(network; producer_growth = g, functional_response) + @test g(1, u, model) == 1 * (1 + f0 * u[2]) * u[1] * (1 - (a_ii * u[1]) / 1) + @test g(2, u, model) == 1 * (1 + f0 * u[1]) * u[2] * (1 - (a_ii * u[2]) / 1) + # Change intra and inter-specific competition. + a_ii, a_ij = rand(2) + g = LogisticGrowth(network; a = (a_ii, a_ij)) + model = ModelParameters(network; producer_growth = g, functional_response) + s1 = a_ii * u[1] + a_ij * u[2] + s2 = a_ii * u[2] + a_ij * u[1] + @test g(1, u, model) == 1 * (1 + f0 * u[2]) * u[1] * (1 - s1 / 1) + @test g(2, u, model) == 1 * (1 + f0 * u[1]) * u[2] * (1 - s2 / 1) + i == 2 && @test g(3, u, model) == 0.0 end - # Test producer competition - # 1 & 2 producer; 3 consumer - A = [0 0 0; 0 0 0; 0 0 1] - foodweb = FoodWeb(A; quiet = true) - rates = BioRates(foodweb; d = 0) - # K = 1, intercompetition (default) - p = ModelParameters( - foodweb; - producer_growth = LogisticGrowth(foodweb; αii = 1.0, αij = 1.0), - biorates = rates - ) - @test simulates(p, [0.5, 0.5, 0.5])[1:2, end] == [0.5, 0.5] + # Test a simple simulation with producer competition. + foodweb = FoodWeb([0 0 0; 0 0 0; 0 0 1]; quiet = true) + biorates = BioRates(foodweb; d = 0) + producer_growth = LogisticGrowth(foodweb; a = (diag = 1, offdiag = 1)) + model = ModelParameters(foodweb; producer_growth, biorates) + u0 = 0.5 # All species have an initial biomass of u0. + sol = simulates(model, u0; verbose = false) + @test sol[1:2, end] == [0.5, 0.5] - p_inter_only = ModelParameters( - foodweb; - producer_growth = LogisticGrowth(foodweb; αii = 0.0, αij = 1.0), - biorates = rates, - ) - s_inter_only = simulates(p_inter_only, [0.5, 0.5, 0.5]) - p_intra_only = ModelParameters( - foodweb; - producer_growth = LogisticGrowth(foodweb; αii = 1.0, αij = 0.0), - biorates = rates, - ) - s_intra_only = simulates(p_intra_only, [0.5, 0.5, 0.5]) - @test s_inter_only[1:2, end] == s_intra_only[1:2, end] ≈ [1.0, 1.0] - - # Extern method with facilitation - foodweb = FoodWeb([0 0 0; 0 0 0; 1 1 0]) # 1 & 2 producers - multiplex_network = MultiplexNetwork(foodweb; C_facilitation = 1.0) - # Facilitation to 0 <=> the growth is unchanged (compared to above section) - multiplex_network.layers[:facilitation].intensity = 0.0 - response = ClassicResponse(multiplex_network) # avoid warning - p = ModelParameters(multiplex_network; functional_response = response) - K, r = p.producer_growth.Kᵢ, p.biorates.r - B = [1, 1, 1] - @test p.producer_growth(1, B, r, multiplex_network) == 0 - @test p.producer_growth(2, B, r, multiplex_network) == 0 - @test p.producer_growth(3, B, r, multiplex_network) == 0 - B = [0.5, 0.5, 0.5] - @test p.producer_growth(1, B, r, multiplex_network) == 0.25 - @test p.producer_growth(2, B, r, multiplex_network) == 0.25 - @test p.producer_growth(3, B, r, multiplex_network) == 0 - rates = BioRates(multiplex_network; r = 2) - p = ModelParameters(multiplex_network; functional_response = response, biorates = rates) - K, r = p.producer_growth.Kᵢ, p.biorates.r - @test p.producer_growth(1, B, r, multiplex_network) == 0.5 - @test p.producer_growth(2, B, r, multiplex_network) == 0.5 - @test p.producer_growth(3, B, r, multiplex_network) == 0 - # Facilitation > 0 <=> the growth is changed (compared to above section) - for f0 in [1.0, 2.0, 5.0, 10.0] - multiplex_network.layers[:facilitation].intensity = f0 - p = ModelParameters(multiplex_network; functional_response = response) - K, r = p.producer_growth.Kᵢ, p.biorates.r - B = [1, 1, 1] - @test p.producer_growth( - 1, - B, - r, - multiplex_network, - ) == 0 - @test p.producer_growth( - 2, - B, - r, - multiplex_network, - ) == 0 - @test p.producer_growth( - 3, - B, - r, - multiplex_network, - ) == 0 - B = [0.5, 0.5, 0.5] - @test p.producer_growth( - 1, - B, - r, - multiplex_network - ) == 0.25 * (1 + f0) - @test p.producer_growth( - 2, - B, - r, - multiplex_network - ) == 0.25 * (1 + f0) - @test p.producer_growth( - 3, - B, - r, - multiplex_network - ) == 0 - rates = BioRates(multiplex_network; r = 2) - p = ModelParameters( - multiplex_network; - functional_response = response, - biorates = rates - ) - K, r = p.producer_growth.Kᵢ, p.biorates.r - @test p.producer_growth( - 1, - B, - r, - multiplex_network, - ) == 0.5 * (1 + f0) - @test p.producer_growth( - 2, - B, - r, - multiplex_network, - ) == 0.5 * (1 + f0) - @test p.producer_growth( - 3, - B, - r, - multiplex_network, - ) == 0 + # Test a simulation with only inter-specific producer competition and + # another simulation with only intra-specific competition. + kwargs = [Dict(:a => (0, 1)), Dict(:a => (1, 0))] + for kw in kwargs + producer_growth = LogisticGrowth(foodweb; kw...) + model = ModelParameters(foodweb; producer_growth, biorates) + sol = simulates(model, u0; verbose = false) + @test sol[1:2, end] ≈ [1, 1] end + end diff --git a/test/model/test-set_temperature.jl b/test/model/test-set_temperature.jl index 2b110c97..7f8a34ea 100644 --- a/test/model/test-set_temperature.jl +++ b/test/model/test-set_temperature.jl @@ -2,26 +2,26 @@ A = [0 0 0; 1 0 0; 0 1 0] foodweb = FoodWeb(A) foodweb.metabolic_class = ["producer", "invertebrate", "ectotherm vertebrate"] foodweb.M = [1.0, 10.0, 10.0] -temp = 303.15 # temperature in Kelvin +temp = 303.15 # Kelvin. @testset "set temperature" begin - p = ModelParameters(foodweb) # bioenergetic default - - # check FR argument error + # By default the functional response is 'bioenergetic', + # activating the temperature dependence should throw an error. + p = ModelParameters(foodweb) @test_throws ArgumentError set_temperature!(p, temp, ExponentialBA()) - # No temperature response + # No temperature response. set_temperature!(p, temp, NoTemperatureResponse()) @test p.environment.T == temp @test p.temperature_response == NoTemperatureResponse() - # MP with Classic Response - p = ModelParameters(foodweb; functional_response = ClassicResponse(foodweb)) - # Exponential Boltzmann-Arrhenius temperature dependence + # Exponential Boltzmann-Arrhenius temperature dependence + p = ModelParameters(foodweb; functional_response = ClassicResponse(foodweb)) set_temperature!(p, temp, ExponentialBA()) @test p.environment.T == temp - @test p.environment.K == exp_ba_vector_rate(foodweb, temp, exp_ba_carrying_capacity()) + @test p.producer_growth.K == + exp_ba_vector_rate(foodweb, temp, exp_ba_carrying_capacity()) @test p.biorates.r == exp_ba_vector_rate(foodweb, temp, exp_ba_growth()) @test p.biorates.x == exp_ba_vector_rate(foodweb, temp, exp_ba_metabolism()) @test p.functional_response.hₜ == @@ -29,7 +29,6 @@ temp = 303.15 # temperature in Kelvin @test p.functional_response.aᵣ == exp_ba_matrix_rate(foodweb, temp, exp_ba_attack_rate()) @test typeof(p.temperature_response) == ExponentialBA - @test p.temperature_response.r == exp_ba_growth() @test p.temperature_response.x == exp_ba_metabolism() @test p.temperature_response.aᵣ == exp_ba_attack_rate() @@ -38,9 +37,9 @@ temp = 303.15 # temperature in Kelvin end @testset "Exponential BA customisation in set_temperature" begin - # MP with Classic Response + + # Exponential Boltzmann-Arrhenius temperature dependence. p = ModelParameters(foodweb; functional_response = ClassicResponse(foodweb)) - # Exponential Boltzmann-Arrhenius temperature dependence set_temperature!( p, temp, diff --git a/test/model/test-simulate.jl b/test/model/test-simulate.jl index cd9d594e..770bc1fd 100644 --- a/test/model/test-simulate.jl +++ b/test/model/test-simulate.jl @@ -1,10 +1,10 @@ @testset "Simulate" begin - # Set up + # Set up. foodweb = FoodWeb([0 0; 1 0]) params = ModelParameters(foodweb; biorates = BioRates(foodweb; d = 0)) - # Solution converges + # Solution converges. solution1 = simulates(params, [0.5, 0.5]; verbose = false) @test is_terminated(solution1) solution2 = simulates(params, [0.3, 0.3]; saveat = 0.25, tmax = 10, verbose = false) @@ -13,16 +13,16 @@ @test is_success(solution3) - # Initial biomass + # Initial biomass. @test solution1.u[begin] == [0.5, 0.5] @test solution2.u[begin] == [0.3, 0.3] @test solution3.u[begin] == [0.2, 0.2] - # Timesteps + # Timesteps. @test all([t ∈ Set(solution2.t) for t in (0:0.25:10)]) @test all([t ∈ Set(solution3.t) for t in (0:0.5:5)]) - # If biomass start at 0, biomass stay at 0 + # If biomass start at 0, biomass stay at 0. solution_null = simulates(params, [0.0, 0.0]; callback = nothing) @test all(hcat(solution_null.u...) .== 0) @test keys(get_extinct_species(solution_null)) == Set([1, 2]) diff --git a/test/model/test-zombies.jl b/test/model/test-zombies.jl index aed6e0aa..edc2847c 100644 --- a/test/model/test-zombies.jl +++ b/test/model/test-zombies.jl @@ -8,7 +8,7 @@ sol = simulates( params, B0; - callback = ExtinctionCallback(1e-5, true), + callback = ExtinctionCallback(1e-5, params, true), tmax = 300_000, verbose = false, compare_rtol = 1e-6, diff --git a/test/runtests.jl b/test/runtests.jl index cbe1ad4a..639b4570 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -42,6 +42,7 @@ test_files = [ "model/test-consumption.jl", "model/test-simulate.jl", "model/test-zombies.jl", + "model/test-nutrient_intake.jl", "model/test-set_temperature.jl", "measures/test-functioning.jl", "measures/test-stability.jl", @@ -54,10 +55,28 @@ test_files = [ function simulates(parms, B0; compare_atol = nothing, compare_rtol = nothing, kwargs...) g = EcologicalNetworksDynamics.simulate(parms, B0; verbose = false, kwargs...) - # Compare with raw specialized code. - xp, data = Logging.with_logger(() -> generate_dbdt(parms, :raw), Logging.NullLogger()) - # Guard against explosive compilation times with this approach. - if SyntaxTree.callcount(xp) <= 20_000 # wild rule of thumb + if is_boostable(parms, :raw) + # Compare with raw specialized code. + xp, data = + Logging.with_logger(() -> generate_dbdt(parms, :raw), Logging.NullLogger()) + # Guard against explosive compilation times with this approach. + if SyntaxTree.callcount(xp) <= 20_000 # wild rule of thumb + dbdt = eval(xp) + s = EcologicalNetworksDynamics.simulate( + parms, + B0; + diff_code_data = (dbdt, data), + verbose = false, + kwargs..., + ) + compare_generic_vs_specialized(g, s, :raw, compare_atol, compare_rtol) + end + end + + if is_boostable(parms, :compact) + # Compare with compact specialized code. + xp, data = + Logging.with_logger(() -> generate_dbdt(parms, :compact), Logging.NullLogger()) dbdt = eval(xp) s = EcologicalNetworksDynamics.simulate( parms, @@ -66,22 +85,9 @@ function simulates(parms, B0; compare_atol = nothing, compare_rtol = nothing, kw verbose = false, kwargs..., ) - compare_generic_vs_specialized(g, s, :raw, compare_atol, compare_rtol) + compare_generic_vs_specialized(g, s, :compact, compare_atol, compare_rtol) end - # Compare with compact specialized code. - xp, data = - Logging.with_logger(() -> generate_dbdt(parms, :compact), Logging.NullLogger()) - dbdt = eval(xp) - s = EcologicalNetworksDynamics.simulate( - parms, - B0; - diff_code_data = (dbdt, data), - verbose = false, - kwargs..., - ) - compare_generic_vs_specialized(g, s, :compact, compare_atol, compare_rtol) - g end diff --git a/use_cases/diversity-vs-nti-intensity.jl b/use_cases/diversity-vs-nti-intensity.jl index 195e67a7..9d047ad8 100644 --- a/use_cases/diversity-vs-nti-intensity.jl +++ b/use_cases/diversity-vs-nti-intensity.jl @@ -25,7 +25,6 @@ n_foodweb = 50 # Replicates of trophic backbones. # `TerminateSteadyState`. tmax = 10_000 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 @@ -62,6 +61,10 @@ Threads.@threads for i in 1:n_foodweb # Parallelize computation if possible. net = MultiplexNetwork(foodweb; [interaction => (L = L_nti, I = intensity)]...) functional_response = ClassicResponse(net; c = i_intra) params = ModelParameters(net; functional_response) + callback = CallbackSet( + ExtinctionCallback(1e-5, params, verbose), + TerminateSteadyState(1e-8, 1e-6), + ) solution = simulate(params, rand(S); tmax, callback) # Save the final diversity at equilibrium. push!(df_thread, [i, interaction, intensity, richness(solution[end])]) @@ -119,3 +122,6 @@ end font = firasans("Medium") # Label font. Label(fig[1:2, 0], "Diversity variation"; font, rotation = pi / 2, width = 0) Label(fig[3, 2:3], "Interaction intensity"; font, height = 0) + +# To save the figure, uncomment and execute the line below. +# save("/tmp/plot.png", fig; resolution = (450, 300), px_per_unit = 3) diff --git a/use_cases/nutrients-competition-exclusion.jl b/use_cases/nutrients-competition-exclusion.jl new file mode 100644 index 00000000..521fd428 --- /dev/null +++ b/use_cases/nutrients-competition-exclusion.jl @@ -0,0 +1,65 @@ +import AlgebraOfGraphics: firasans, set_aog_theme! +using CairoMakie +using EcologicalNetworksDynamics +using Random + +foodweb = FoodWeb([0 0; 0 0]) # Two plants competing for nutrients. +half_saturation = [0.3 0.9; 0.9 0.3] # Row <-> plants & column <-> nutrients. +turnover = 0.9 +supply = 10 +concentration = 1 +d = [0.6, 1.2] +r = [1, 2] +tmax = 50 +extinction_threshold = 1e-5 + +sol_vec = [] +for n_nutrients in 1:2 + producer_growth = NutrientIntake( + foodweb; + n_nutrients, + supply, + half_saturation = half_saturation[:, 1:n_nutrients], + turnover, + concentration, + ) + biorates = BioRates(foodweb; d, r) + p = ModelParameters(foodweb; producer_growth, biorates) + Random.seed!(113) # Set seed for reproducibility of initial conditions. + N0 = 0.1 .+ 3 * rand(n_nutrients) + B0 = [0.5, 0.5] + callback = ExtinctionCallback(extinction_threshold, p, true) + sol = simulate(p, B0; N0, tmax, alg_hints = [:stiff], reltol = 1e-5, callback) + push!(sol_vec, sol) +end + +set_aog_theme!() # AlgebraOfGraphics theme. +fig = Figure() + +ax1 = Axis(fig[1, 1]; xlabel = "", ylabel = "", title = "1 nutrient: exclusion") +sol = sol_vec[1] +plant1_line = lines!(sol.t, sol[1, :]; color = :red) +plant2_line = lines!(sol.t, sol[2, :]; color = :green) +nutrient1_line = lines!(sol.t, sol[3, :]; color = :blue, linestyle = :dot) + +ax2 = Axis(fig[1, 2]; xlabel = "", ylabel = "", title = "2 nutrients: coexistence") +sol = sol_vec[2] +plant1_line = lines!(sol.t, sol[1, :]; color = :red) +plant2_line = lines!(sol.t, sol[2, :]; color = :green) +nutrient1_line = lines!(sol.t, sol[3, :]; color = :blue, linestyle = :dot) +nutrient2_line = lines!(sol.t, sol[4, :]; color = :black, linestyle = :dot) + +font = firasans("Medium") # Label font. +Label(fig[1, 0], "Biomass"; font, rotation = pi / 2, width = 0, tellheight = false) +Label(fig[2, 1:2], "Time"; font, height = 0, tellheight = true) +linkyaxes!(ax1, ax2) +hideydecorations!(ax2; ticks = false) +Legend( + fig[0, :], + [plant1_line, plant2_line, nutrient1_line, nutrient2_line], + ["Plant 1", "Plant 2", "Nutrient 1", "Nutrient 2"]; + orientation = :horizontal, + tellheight = true, # Adjust the height of the legend sub-figure. + tellwidth = false, # Do not adjust width of the orbit diagram. +) +save("/tmp/plot.png", fig; resolution = (450, 300), px_per_unit = 3) diff --git a/use_cases/paradox-enrichment.jl b/use_cases/paradox-enrichment.jl index 0ca6f065..bdd1e5e1 100644 --- a/use_cases/paradox-enrichment.jl +++ b/use_cases/paradox-enrichment.jl @@ -29,8 +29,8 @@ df = DataFrame(; # Run simulations: compute equilibirum biomass for each carrying capacity. @info "Start simulations..." for K in K_values - environment = Environment(foodweb; K) - params = ModelParameters(foodweb; functional_response, environment) + producer_growth = LogisticGrowth(foodweb; K) + params = ModelParameters(foodweb; functional_response, producer_growth) B0 = rand(2) # Inital biomass. solution = simulate(params, B0; tmax, verbose) extrema = biomass_extrema(solution, "10%") @@ -61,3 +61,6 @@ Legend( tellheight = true, # Adjust the height of the legend sub-figure. tellwidth = false, # Do not adjust width of the orbit diagram. ) + +# To save the figure, uncomment and execute the line below. +# save("/tmp/plot.png", fig; resolution = (450, 300), px_per_unit = 3) diff --git a/use_cases/paradox_enrichment_nutrientintake.jl b/use_cases/paradox_enrichment_nutrientintake.jl index 7dea5c8d..e3e6c540 100644 --- a/use_cases/paradox_enrichment_nutrientintake.jl +++ b/use_cases/paradox_enrichment_nutrientintake.jl @@ -17,8 +17,8 @@ end foodweb = FoodWeb([2 => 1]); # 2 eats 1. functional_response = ClassicResponse(foodweb; aᵣ = 1, hₜ = 1, h = 1); -S_values = LinRange(1, 60, 60) -tmax = 1_000 # Simulation length. +S_values = LinRange(0, 70, 20) +tmax = 10_000 # Simulation length. verbose = false # Do not show '@info' messages during the simulation. df = DataFrame(; S = Float64[], @@ -30,23 +30,31 @@ df = DataFrame(; # Run simulations: compute equilibirum biomass for each carrying capacity. @info "Start simulations..." -for s in S_values - for i in 1:20 - k1 = round(rand(Uniform(0.1, 0.2), 1)[1], digits = 2) - k2 = round(rand(Uniform(0.1, 0.2), 1)[1], digits = 2) - d = round(rand(Uniform(0.1, 0.4), 1)[1], digits = 2) - growthmodel = NutrientIntake(foodweb, - supply = s, - half_saturation = hcat(k1,k2), - turnover = d +for supply in S_values + for _ in 1:10 + half_saturation = 0.1 + turnover = 0.1 + producer_growth = NutrientIntake(foodweb; supply, half_saturation, turnover) + params = ModelParameters(foodweb; functional_response, producer_growth) + callback = ExtinctionCallback(1e-6, params, verbose) + B0 = 1 .+ 3 * rand(2) # Inital biomass. + N0 = 1 .+ 3 * rand(2) # Initial nutrient abundances. + solution = simulate( + params, + B0; + N0, + tmax, + verbose, + alg_hints = [:stiff], + callback, + reltol = 1e-5, + ) + if solution.retcode == ReturnCode.Success + extrema = biomass_extrema(solution, "10%") + push!( + df, + [supply, extrema[1].min, extrema[1].max, extrema[2].min, extrema[2].max], ) - params = ModelParameters(foodweb; functional_response, producer_growth = growthmodel) - B0 = rand(2) # Inital biomass. - N0 = rand(2) - solution = simulate(params, B0; N0 = N0, tmax, verbose, alg_hints=[:stiff]) - extrema = biomass_extrema(solution, "10%") - if solution.retcode == ReturnCode.Terminated - push!(df, [s, extrema[1].min, extrema[1].max, extrema[2].min, extrema[2].max]) end end @info "Simulation for supply S = $s done." @@ -55,20 +63,24 @@ end # Plot the orbit diagram with Makie. df2 = groupby(df, :S) -df2 = combine(df2, - :B_resource_min => mean, - :B_resource_max => mean, - :B_consumer_min => mean, - :B_consumer_max => mean) +df2 = combine( + df2, + :B_resource_min => mean, + :B_resource_max => mean, + :B_consumer_min => mean, + :B_consumer_max => mean, +) 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 = "Supply, S", ylabel = "Equilibrium biomass") -resource_line = scatterlines!(df2.S, df2.B_resource_min_mean; color = c_r, markercolor = c_r) +resource_line = + scatterlines!(df2.S, df2.B_resource_min_mean; color = c_r, markercolor = c_r) scatterlines!(df2.S, df2.B_resource_max_mean; color = c_r, markercolor = c_r) -consumer_line = scatterlines!(df2.S, df2.B_consumer_min_mean; color = c_c, markercolor = c_c) +consumer_line = + scatterlines!(df2.S, df2.B_consumer_min_mean; color = c_c, markercolor = c_c) scatterlines!(df2.S, df2.B_consumer_max_mean; color = c_c, markercolor = c_c) Legend( fig[1, 1], @@ -78,4 +90,6 @@ Legend( tellheight = true, # Adjust the height of the legend sub-figure. tellwidth = false, # Do not adjust width of the orbit diagram. ) -fig \ No newline at end of file + +# To save the figure, uncomment and execute the line below. +# save("/tmp/plot.png", fig; resolution = (450, 300), px_per_unit = 3) diff --git a/use_cases/persistence-vs-producer-competition.jl b/use_cases/persistence-vs-producer-competition.jl index 03c14c3f..056b2e44 100644 --- a/use_cases/persistence-vs-producer-competition.jl +++ b/use_cases/persistence-vs-producer-competition.jl @@ -31,14 +31,18 @@ end # Main simulation loop. # 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 (i_C, C) in C_values # Parallelize on connctance values. +Threads.@threads for i_C in 1:length(C_values) # Parallelize on connctance values. + C = C_values[i_C] 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 - producer_competition = ProducerCompetition(foodweb; αii, αij) - environment = Environment(foodweb; K = standardize_K(foodweb, K, αij)) - params = ModelParameters(foodweb; producer_competition, environment) + producer_growth = LogisticGrowth( + foodweb; + K = standardize_K(foodweb, K, αij), + a = (diag = αii, offdiag = αij), + ) + params = ModelParameters(foodweb; producer_growth) B0 = rand(S) # Initial biomass. solution = simulate(params, B0; extinction_threshold, tmax, verbose) # Measure species persistence i.e. the number of species that have @@ -90,3 +94,6 @@ Legend( tellheight = true, # Adjust top subfigure height to legend height. tellwidth = false, # Do not adjust bottom subfigure width to legend width. ) + +# To save the figure, uncomment and execute the line below. +# save("/tmp/plot.png", fig; resolution = (450, 300), px_per_unit = 3) diff --git a/use_cases/tritrophic-biomass-vs-temperature.jl b/use_cases/tritrophic-biomass-vs-temperature.jl index b2c8a877..58530924 100644 --- a/use_cases/tritrophic-biomass-vs-temperature.jl +++ b/use_cases/tritrophic-biomass-vs-temperature.jl @@ -25,14 +25,14 @@ df = DataFrame(; T = Float64[], trophic_level = Int64[], Beq = Float64[]) # Simulation length and time steps are taken from Binzer et al., 2012. tmax = 315_360_000_000 # From Binzer et al., 2012. verbose = false # Do not show '@info' messages during simulation run. +callback = ExtinctionCallback(1e-6, params, verbose) # Remove TerminateSteadyState callback. # Run simulations for each temperature across gradient. @info "Start simulations..." 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. - callback = ExtinctionCallback(1e-6, verbose) # Remove TerminateSteadyState callback. + B0 = params.producer_growth.K[1] / 8 # Inital biomass. solution = simulate(params, [B0]; tmax, callback) Beq_vec = solution[end] for (Beq, tlvl) in zip(Beq_vec, trophic_lvl) @@ -63,3 +63,6 @@ Legend( tellheight = true, # Adjust top subfigure height to legend height. tellwidth = false, # Do not adjust bottom subfigure width to legend width. ) + +# To save the figure, uncomment and execute the line below. +# save("/tmp/plot.png", fig; resolution = (450, 300), px_per_unit = 3) diff --git a/use_cases/trophic-classes-vs-mortality.jl b/use_cases/trophic-classes-vs-mortality.jl index aadea3d3..d774e415 100644 --- a/use_cases/trophic-classes-vs-mortality.jl +++ b/use_cases/trophic-classes-vs-mortality.jl @@ -27,7 +27,6 @@ check_cycle = true # Check that generated food webs do not contain cycle(s). # To do so, we have lowered the `abstol` and `reltol` arguments of # `TerminateSteadyState`. verbose = false # Do not show '@info' messages during simulation run. -callback = CallbackSet(ExtinctionCallback(1e-5, verbose), TerminateSteadyState(1e-8, 1e-6)) tmax = 10_000 # Main simulation loop. @@ -55,10 +54,12 @@ Threads.@threads for i in 1:n_foodweb # Parallelize computation when possible. d = d0 .* allometric_rate(foodweb, DefaultMortalityParams()) biorates = BioRates(foodweb; d) p = ModelParameters(foodweb; functional_response, biorates) + callback = CallbackSet( + ExtinctionCallback(1e-5, p, verbose), + TerminateSteadyState(1e-8, 1e-6), + ) # Prepare simulation boost. - dBdt_expr, data = generate_dbdt(p, :compact) - dBdt! = eval(dBdt_expr) - solution = simulate(p, B0; tmax, callback, diff_code_data = (dBdt!, data)) + solution = simulate(p, B0; tmax, callback) extinct_sp = keys(get_extinct_species(solution)) surviving_sp = setdiff(1:richness(foodweb), extinct_sp) # Species classes are given by the the dynamics with zero mortality rate. @@ -135,3 +136,6 @@ Legend( tellheight = true, # Adjust top subfigure height to legend height. tellwidth = false, # Do not adjust bottom subfigure width to legend width. ) + +# To save the figure, uncomment and execute the line below. +# save("/tmp/plot.png", fig; resolution = (450, 300), px_per_unit = 3)