Rewrote insect submodel

It is now contained within a single function that can be called at
need, rather than being calculated every day for every pixel.

src/core/simulation.jl

@@ -118,7 +118,6 @@ Execute one update of the model.
 function stepsimulation!(model::AgentBasedModel)
     with_logger(model.logger) do
         @info "Simulating day $(model.date)."
-        updateinsectbiomass!(model)
         for a in Schedulers.ByType((Farmer,FarmPlot,Animal), true)(model)
             try #The animal may have been killed
                 stepagent!(model[a], model)

src/nature/insects.jl

 ### Persefone - a socio-economic-ecological model of European agricultural landscapes.
 ###
-### This file contains the functions for calculating insect biomass.
+### This file contains submodel that calculates insect biomass
 ###
 
-##FIXME Updating every pixel every day is really expensive! I should seriously
-## think about switching this to a "calculate-on-demand" design.
-## (The problem is that that works for all models except the pesticide one, which
-## is a Markov Chain model that cannot be calculated without knowing the previous
-## value at the location.)
-
 ##TODO write tests for this
 
 """
-    updateinsectbiomass!(model)
-
-Calculate a new insect biomass value for each pixel in the landscape.
-"""
-function updateinsectbiomass!(model::AgentBasedModel)
-    @debug "Updating insect biomass"
-    if @param(nature.insectmodel) == "seasonal"
-        updatefunction = seasonalmodel
-    elseif @param(nature.insectmodel) == "seasonalhabitat"
-        updatefunction = seasonalhabitatmodel
-    else
-        Base.error("Unknown insect model $(@param(nature.insectmodel)).")
-        #TODO replace with exception
-    end
-    for p in model.landscape
-        p.insectbiomass = updatefunction(p, model)
-    end
-end
+    insectbiomass(pixel, model)
 
+Calculate the insect biomass in this location, using the factors configured
+in the `nature.insectmodel` settings (any of: "season", "habitat", "weather",
+"pesticides"). Returns a float value in mg/m².
 
-"""
-    seasonalmodel(pixel)
+**Biological note:** this is a very approximate calculation! Insect biomass
+varies wildly in time and space and is hard to measure. This calculation is
+based on the idea of a parabolic seasonal development of insect abundance,
+modified habitat suitability, weather, and pesticide application. Although it
+is based on empirical studies, it can only deliver a rough, order-of-magnitude
+estimation of likely insect biomass in a given location.
 
-Calculate an insect biomass value based purely on the calendar day, following
-a parabolic function as per Paquette et al. (2013). Assumes a maximum insect
-density of 300 mg/m² (cf. Odderskaer et al. 1997, Püttmanns et al. 2022).
+**Sources:**
 
-*[Note that this is a very approximate calculation!]*
+- Odderskær et al. (1997). Skylark Reproduction in Pesticide Treated and
+Untreated Fields (32; Pesticides Research). Danish Environmental Protection Agency.
+- Grüebler et al. (2008). A predictive model of the density of airborne
+insects in agricultural environments. Agriculture, Ecosystems & Environment,
+123(1), 75–80. https://doi.org/10.1016/j.agee.2007.05.001
+- Paquette et al. (2013). Seasonal patterns in Tree Swallow prey (Diptera)
+abundance are affected by agricultural intensification. Ecological Applications,
+23(1), 122–133. https://doi.org/10.1890/12-0068.1
+- Püttmanns et al. (2022). Habitat use and foraging parameters of breeding
+Skylarks indicate no seasonal decrease in food availability in heterogeneous
+farmland. Ecology and Evolution, 12(1), e8461. https://doi.org/10.1002/ece3.8461
 """
-function seasonalmodel(pixel::Pixel, model::AgentBasedModel)
-    calendarday = dayofyear(model.date)
-    biomass = -0.085(calendarday-187)^2+300
-    biomass > 0 ? biomass : 0
-end
+function insectbiomass(pixel::Pixel, model::AgentBasedModel)::Float64
 
-"""
-    seasonalhabitatmodel(pixel)
+    ## if no factors are configure, insect abundance defaults to 300 mg/m²,
+    ## a value in the upper range of insect biomass density in agricultural landscapes
+    baseline = 300.0
+    seasonfactor = 0.0
+    habitatfactor = 1.0
+    weatherfactor = 1.0
+    pesticidefactor = 1.0
 
-Calculate an insect biomass value based on the calendar day and habitat type,
-following a parabolic function.
+    ## parabolic curve of seasonal development,
+    ## based on fig. 3a in Paquette et al. (2013)
+    if "season" in @param(nature.insectmodel)
+        calendarday = dayofyear(model.date)
+        seasonfactor = -0.085(calendarday-187)^2
+    end
 
-See Odderskaer et al. (1997), Grübler et al. (2008), Paquette et al. (2013),
-and Püttmanns et al. (2022).
+    ## habitat dependence of insect biomass,
+    ## based on fig. 1 in Grübler et al. (2008) and fig. 5c in Püttmanns et al. (2022)
+    if "habitat" in @param(nature.insectmodel)
+        if pixel.landcover == soil || pixel.landcover == forest
+            habitatfactor = 2.0
+        elseif pixel.landcover == grass
+            habitatfactor = 1.5
+        elseif pixel.landcover == agriculture
+            habitatfactor = 1.0
+        else
+            habitatfactor = 0.0
+        end
+    end
 
-*[Note that this is a very approximate calculation!]*
-"""
-function seasonalhabitatmodel(pixel::Pixel, model::AgentBasedModel)
-    calendarday = dayofyear(model.date)
-    habitatfactor = 0.0
-    if pixel.landcover == soil || pixel.landcover == forest
-        habitatfactor = 2.0
-    elseif pixel.landcover == grass
-        habitatfactor = 1.5
-    elseif pixel.landcover == agriculture
-        habitatfactor = 1.0
+    ## temperature dependence of insect biomass,
+    ## based on fig. 3b in Paquette et al. (2013)
+    ##XXX (and possibly table 3 in Grübler et al. (2008))
+    if "weather" in @param(nature.insectmodel)
+        #TODO add this once we have a working weather module
+    end
+
+    ## effect of pesticides on insect abundance,
+    ## based on figs. 3.6 and 3.7 in Odderskær et al. (1997)
+    ## Note that this is a simplification: it ignores that insect biomass
+    ## tends to rise rapidly a few weeks after pesticide application.
+    if "pesticides" in @param(nature.insectmodel) && pesticide in pixel.events
+        pesticidefactor = 0.5
     end
-    biomass = -0.085(calendarday-187)^2+300*habitatfactor
-    biomass > 0 ? biomass : 0
-end
 
-#TODO add models that include weather and pesticide application
+    ## calculate biomass using a parabolic equation in the vertex form
+    biomass = seasonfactor+baseline*habitatfactor*pesticidefactor*weatherfactor
+    biomass > 0 ? biomass : 0.0
+end

src/parameters.toml

@@ -29,7 +29,7 @@ farmmodel = "FieldManager" # which version of the farm model to use (not yet imp
 targetspecies = ["Wolpertinger", "Wyvern"] # list of target species to simulate
 popoutfreq = "daily" # output frequency population-level data, daily/monthly/yearly/end/never
 indoutfreq = "end" # output frequency individual-level data, daily/monthly/yearly/end/never
-insectmodel = "seasonalhabitat" # how to calculate insect biomass, seasonal/seasonalhabitat/full
+insectmodel = ["season", "habitat", "pesticides"] # which factors affect insect growth ("weather" is not yet implemented)
 	
 [crop]
 cropmodel = "linear" # crop growth model to use, "linear" or "aquacrop" (not yet implemented)

src/world/landscape.jl

@@ -21,7 +21,6 @@ mutable struct Pixel
     landcover::LandCover
     fieldid::Union{Missing,Int64}
     events::Vector{EventType}
-    insectbiomass::Float64 # in mg/m²
 end
 
 """
@@ -58,7 +57,7 @@ function initlandscape(landcovermap::String, farmfieldsmap::String)
             lcv = LandCover(Int(lc/10))
             ff = farmfields[x,y][1]
             !(ismissing(ff)) && (ff = Int64(ff))
-            landscape[x,y] = Pixel(lcv, ff, Vector{Symbol}(), 0.0)
+            landscape[x,y] = Pixel(lcv, ff, Vector{Symbol}())
         end
     end
     return landscape

test/nature_tests.jl

@@ -6,7 +6,7 @@
 @testset "Habitat macros" begin
     # set up the testing landscape
     model = smalltestlandscape()
-    model.landscape[6,6] = Pixel(Ps.agriculture, 1, [], 0.0)
+    model.landscape[6,6] = Pixel(Ps.agriculture, 1, [])
     species::Dict{String,Any} = Dict("name"=>"test_animal")
     add_agent!((6,6), FarmPlot, model, [(6,6)], Ps.wheat, 1.2, 3.4)
     add_agent!((3,3), Animal, model, species, Ps.male, 1)
@@ -136,3 +136,11 @@ end
     @test typeof(@rand()) == Float64
     @test @rand([true, true])
 end
+
+
+@testset "Insect submodel" begin
+    model = smalltestlandscape()
+
+    #TODO create a set of pixels and check if their calculated biomass
+    # meets the expected values
+end

test/runtests.jl

@@ -37,10 +37,10 @@ function smalltestlandscape()
     for x in 1:6
         for y in 1:6
             (x in (1:2) || y in (1:4)) ? lc = Ps.forest : lc = Ps.grass
-            landscape[x,y] = Pixel(lc, missing, [], 0.0)
+            landscape[x,y] = Pixel(lc, missing, [])
         end
     end
-    landscape[6,4] = Pixel(Ps.water, 0, [], 0.0)
+    landscape[6,4] = Pixel(Ps.water, 0, [])
     space = GridSpace(size(landscape), periodic=false)
     properties = Dict{Symbol,Any}(:date=>TESTSETTINGS["core.startdate"],
                                   :landscape=>landscape,