Incorporate Iago's suggestions

src/BEFWM2.jl

@@ -48,6 +48,7 @@ export A_interference_full
 export A_refuge_full
 export allometric_rate
 export AllometricParams
+export attack_rate
 export BioEnergeticFunctionalResponse
 export BioenergeticResponse
 export BioRates
@@ -63,6 +64,7 @@ export DefaultMortalityParams
 export draw_asymmetric_links
 export draw_symmetric_links
 export Environment
+export efficiency
 export find_steady_state
 export fitin
 export FoodWeb
@@ -74,6 +76,7 @@ export foodweb_simpson
 export FunctionalResponse
 export FunctionalResponse
 export generate_dbdt
+export handling_time
 export homogeneous_preference
 export interaction_names
 export Layer

src/inputs/biological_rates.jl

@@ -193,16 +193,17 @@ function BioRates(
     e = efficiency(network),
 )
     S = richness(network)
-    Rates = [d, r, x, y]
+    rate_list = [d, r, x, y]
 
-    # Perform sanity checks and vectorize rates if necessary
-    for (i, rate) in enumerate(Rates)
-        isa(rate, Real) ? (Rates[i] = fill(rate, S)) : @check_equal_richness length(rate) S
+    # Perform sanity checks and vectorize rate if necessary
+    for (i, rate) in enumerate(rate_list)
+        isa(rate, Real) ? (rate_list[i] = fill(rate, S)) :
+        @check_equal_richness length(rate) S
     end
     @check_size_is_richness² e S
 
     # Output
-    d, r, x, y = Rates
+    d, r, x, y = rate_list
     BioRates(d, r, x, y, e)
 end
 

src/inputs/functional_response.jl

@@ -93,7 +93,7 @@ The efficiency depends on the metabolic class of the prey:
 - if prey is producter, efficiency is `e_herbivore`
 - otherwise efficiency is `e_carnivore`
 
-Default values are taken from *add ref*.
+Default values are taken from Miele et al. 2019 (PLOS Comp.).
 """
 function efficiency(net::EcologicalNetwork; e_herb = 0.45, e_carn = 0.85)
     S = richness(net)
@@ -266,7 +266,8 @@ and [`FunctionalResponse`](@ref).
 function (F::ClassicResponse)(B, i, j, mᵢ)
     num = F.ω[i, j] * F.aᵣ[i, j] * B[j]^F.h
     denom = 1 + (F.c[i] * B[i]) + sum(F.aᵣ[i, :] .* F.hₜ[i, :] .* F.ω[i, :] .* (B .^ F.h))
-    num / (mᵢ * denom)
+    denom *= mᵢ
+    num / denom
 end
 # Code generation version (raw) (↑ ↑ ↑ DUPLICATED FROM ABOVE ↑ ↑ ↑).
 # (update together as long as the two coexist)
@@ -358,8 +359,9 @@ function (F::ClassicResponse)(B, i, j, aᵣ, network::MultiplexNetwork)
     i0 = network.layers[:interference].intensity
     predator_interfering = A_interference[:, i]
     denom += i0 * sum(B .* predator_interfering)
+    denom *= network.M[i]
 
-    num / (network.M[i] * denom)
+    num / denom
 end
 
 """

src/model/simulate.jl

@@ -58,9 +58,9 @@ 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> biorates = BioRates(foodweb, d=0); # set natural death rate to 0
+julia> br = BioRates(foodweb, d=0); # set natural death rate to 0
 
-julia> params = ModelParameters(foodweb, biorates=biorates); 
+julia> params = ModelParameters(foodweb, biorates=br); 
 
 julia> B0 = [0.5, 0.5]; # set initial biomass
 

test/inputs/test-functional_response.jl

@@ -1,12 +1,12 @@
 @testset "Assimilation efficiency" begin
     foodweb = FoodWeb([0 0 0; 1 0 0; 1 1 0])
     e_expect = sparse([0 0 0; 1 0 0; 1 2 0])
-    @test BEFWM2.efficiency(foodweb; e_herb = 1, e_carn = 2) == e_expect
+    @test efficiency(foodweb; e_herb = 1, e_carn = 2) == e_expect
     e_expect = sparse([0 0 0; 3 0 0; 3 4 0])
-    @test BEFWM2.efficiency(foodweb; e_herb = 3, e_carn = 4) == e_expect
+    @test efficiency(foodweb; e_herb = 3, e_carn = 4) == e_expect
     foodweb = FoodWeb([0 0 0; 1 1 0; 1 1 1])
     e_expect = sparse([0 0 0; 1 2 0; 1 2 2])
-    @test BEFWM2.efficiency(foodweb; e_herb = 1, e_carn = 2) == e_expect
+    @test efficiency(foodweb; e_herb = 1, e_carn = 2) == e_expect
 end
 
 A_2sp = [0 0; 1 0] # 2 eats 1
@@ -60,8 +60,8 @@ end
     foodweb_default = FoodWeb([0 0 0; 1 0 0; 0 1 0]; Z = 10)
     Fclassic_1 = ClassicResponse(foodweb_default)
     @test Fclassic_1.h == 2.0
-    @test Fclassic_1.aᵣ == BEFWM2.attack_rate(foodweb_default)
-    @test Fclassic_1.hₜ == BEFWM2.handling_time(foodweb_default)
+    @test Fclassic_1.aᵣ == attack_rate(foodweb_default)
+    @test Fclassic_1.hₜ == handling_time(foodweb_default)
     @test Fclassic_1.c == [0.0, 0.0, 0.0]
 
     # Custom
@@ -314,25 +314,29 @@ end
 @testset "Generation of default feeding rates" begin
     # All body masses set to 1, expect 0.3 for each trophic interaction
     foodweb = FoodWeb([0 0 0; 1 0 0; 0 1 0]; M = [1, 1, 1])
-    @test BEFWM2.handling_time(foodweb) == [0 0 0; 0.3 0 0; 0 0.3 0]
-    @test BEFWM2.handling_time(foodweb |> MultiplexNetwork) == [0 0 0; 0.3 0 0; 0 0.3 0]
-    @test BEFWM2.attack_rate(foodweb) == [0 0 0; 50 0 0; 0 50 0]
-    @test BEFWM2.attack_rate(foodweb |> MultiplexNetwork) == [0 0 0; 50 0 0; 0 50 0]
+    @test handling_time(foodweb) == [0 0 0; 0.3 0 0; 0 0.3 0]
+    @test handling_time(foodweb |> MultiplexNetwork) == [0 0 0; 0.3 0 0; 0 0.3 0]
+    @test attack_rate(foodweb) == [0 0 0; 50 0 0; 0 50 0]
+    @test attack_rate(foodweb |> MultiplexNetwork) == [0 0 0; 50 0 0; 0 50 0]
 
     # Different body masses, expect different values
     foodweb = FoodWeb([0 0 0; 1 0 0; 0 1 0]; M = [2, 10, 100])
+    a = 0.3 * 10^(-0.48) * 2^(-0.66) # expected ht[2,1]
+    b = 0.3 * 100^(-0.48) * 10^(-0.66) # expected ht[3,2]
     ht_expected = [
         0 0 0
-        0.3*10^(-0.48)*2^(-0.66) 0 0
-        0 0.3*100^(-0.48)*10^(-0.66) 0
+        a 0 0
+        0 b 0
     ]
+    c = 50 * 10^(0.45) # expected ar[2,1]
+    d = 50 * 100^(0.45) * 10^(0.15) # expected ar[3,2]
     ar_expected = [
         0 0 0
-        50*10^(0.45) 0 0
-        0 50*100^(0.45)*10^(0.15) 0
+        c 0 0
+        0 d 0
     ]
-    @test BEFWM2.handling_time(foodweb) ≈ ht_expected
-    @test BEFWM2.handling_time(foodweb |> MultiplexNetwork) ≈ ht_expected
-    @test BEFWM2.attack_rate(foodweb) ≈ ar_expected
-    @test BEFWM2.attack_rate(foodweb |> MultiplexNetwork) ≈ ar_expected
+    @test handling_time(foodweb) ≈ ht_expected
+    @test handling_time(foodweb |> MultiplexNetwork) ≈ ht_expected
+    @test attack_rate(foodweb) ≈ ar_expected
+    @test attack_rate(foodweb |> MultiplexNetwork) ≈ ar_expected
 end