Merge branch 'next' of https://github.com/PoisotLab/BioEnergeticFoodW…

…ebs.jl into next

src/dBdt.jl

@@ -147,6 +147,12 @@ function consumption(parameters, biomass)
 
 end
 
+function density_dependent_mortality(parameters, biomass)
+  mortality_c = parameters[:dc](biomass) .* Int.(.!parameters[:is_producer])
+  mortality_p = parameters[:dp](biomass) .* Int.(parameters[:is_producer])
+  mortality = mortality_c .+ mortality_p
+end
+
 """
 **Derivatives**
 
@@ -172,10 +178,13 @@ function dBdt(derivative, biomass, parameters::Dict{Symbol,Any}, t)
   # Growth
   growth, G = BioEnergeticFoodWebs.get_growth(parameters, biomass; c = nutrients)
 
+  # Mortality
+  mortality = BioEnergeticFoodWebs.density_dependent_mortality(parameters, biomass)
+
   # Balance
   dbdt = zeros(eltype(biomass), length(biomass))
   for i in eachindex(dbdt)
-    dbdt[i] = growth[i] + gain[i] - loss[i]
+    dbdt[i] = growth[i] + gain[i] - loss[i] - mortality[i]
   end
 
   parameters[:productivity] == :nutrients && append!(dbdt, BioEnergeticFoodWebs.nutrientuptake(parameters, biomass, nutrients, G))

src/make_parameters.jl

@@ -26,6 +26,8 @@ matrix A. See documentation for more information. Specifically, the default valu
 | scale_metabolism  | false         | whether to normalize metabolic rates by the growth rate of the smallest producer            |
 | scale_maxcons     | false         | whether to normalize max. consumption rates by metabolic rates                              |
 | productivity      | :species      | type of productivity regulation                                                             |
+| dc                | x -> x .* 0.0 | density dependent mortality function for consumers                                          |
+| dp                | x -> x .* 0.0 | density dependent mortality function for producers                                          |
 | rewire_method     | :none         | method for rewiring the foodweb following extinction events                                 |
 | adbm_trigger      | :extinction   | (ADBM) trigger for ADBM rewiring (on extinctions or periodic with :interval)                |
 | adbm_interval     | 100           | (ADBM) Δt for periodic rewiring                                                             |
@@ -107,6 +109,8 @@ function model_parameters(A;
         bodymass::Array{Float64, 1}=[0.0],
         scale_bodymass::Bool=true,
         vertebrates::Array{Bool, 1}=[false],
+        dc::Function= (x -> x .* 0.0),
+        dp::Function= (x -> x .* 0.0),
         rewire_method::Symbol = :none,
         adbm_trigger::Symbol = :extinction,
         adbm_interval::Int64 = 100,
@@ -321,11 +325,14 @@ function model_parameters(A;
   # Step 18 -- Efficiency matrix
   get_efficiency(parameters)
 
-  # Final Step -- store the parameters in the dict. p
   parameters[:Γh] = parameters[:Γ] .^ parameters[:h]
   parameters[:np] = sum(parameters[:is_producer])
   parameters[:ar] = attack_r
 
+  # Step  19 -- Density dependent mortality
+  parameters[:dc] = dc
+  parameters[:dp] = dp
+
   check_parameters(parameters)
 
   return parameters

test/mortality.jl

@@ -0,0 +1,47 @@
+module TestDDMortality
+    using BioEnergeticFoodWebs
+    using Test
+
+    A = [0 1 0 ; 0 0 1 ; 0 0 0]
+    ddm = (x -> x .* 0.2)
+    p = model_parameters(A, dc = ddm)
+
+    b = [0.5, 0.6, 0.8]
+    PI_death = BioEnergeticFoodWebs.density_dependent_mortality(p, b)
+    @test all(PI_death .== ddm(b) .* Int.(.!p[:is_producer]))
+
+    dbdt = BioEnergeticFoodWebs.dBdt(zeros(3), b, p, 1)
+    # Consumption
+    gain, loss = BioEnergeticFoodWebs.consumption(p, b)
+    # Growth
+    growth, G = BioEnergeticFoodWebs.get_growth(p, b)
+    # Balance
+    balance = zeros(eltype(b), length(b))
+    for i in eachindex(balance)
+      balance[i] = growth[i] + gain[i] - loss[i] - PI_death[i]
+    end
+
+    @test all(dbdt .== balance)
+
+    ddm_prod = (x -> x .^ 2 .* 0.1)
+    p = model_parameters(A, dc = ddm, dp = ddm_prod)
+    PI_death = BioEnergeticFoodWebs.density_dependent_mortality(p, b)
+    mc = ddm(b).* Int.(.!p[:is_producer])
+    mp = ddm_prod(b).* Int.(p[:is_producer])
+    @test all(PI_death .== mc .+ mp)
+
+    dbdt = BioEnergeticFoodWebs.dBdt(zeros(3), b, p, 1)
+    # Consumption
+    gain, loss = BioEnergeticFoodWebs.consumption(p, b)
+    # Growth
+    growth, G = BioEnergeticFoodWebs.get_growth(p, b)
+    # Balance
+    balance = zeros(eltype(b), length(b))
+    for i in eachindex(balance)
+      balance[i] = growth[i] + gain[i] - loss[i] - PI_death[i]
+    end
+
+    @test all(dbdt .== balance)
+
+
+end

test/runtests.jl

@@ -17,7 +17,8 @@ test_files = [
     "biological_rates.jl",
     "temperature_size.jl",
     "productivity.jl",
-    "extinctions.jl"
+    "extinctions.jl",
+    "mortality.jl"
 ]
 
 test_n = 1