From 4e876a092214445f09f32929705f5d745a7834d1 Mon Sep 17 00:00:00 2001
From: Daniel Vedder <daniel.vedder@idiv.de>
Date: Thu, 18 Jul 2024 16:03:44 +0200
Subject: [PATCH] Wrote skylark territory searching functions

---
 CHANGELOG.md                  |   6 +-
 src/crop/farmplot.jl          |   4 +-
 src/nature/individuals.jl     |  13 +--
 src/nature/macros.jl          |   1 -
 src/nature/nature.jl          |   2 +-
 src/nature/populations.jl     |   4 +-
 src/nature/species/skylark.jl | 149 +++++++++++++++++++++++++---------
 7 files changed, 129 insertions(+), 50 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index aa0f220..2f9f15e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,10 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 
+## [1.0.0] - in planning
+
+*Aim: 3 species, 2 crop growth models, farm model, GAEC scenarios, experimental analysis*
+
 ## [0.6.0] - in planning
 
 *Plan: integrate multiple crop growth models, set up first experiments*
@@ -33,7 +37,6 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - `EventType` renamed to `Management` for clarity
 
-- `insectbiomass()` uses units
 
 ### Deprecated
 
@@ -45,6 +48,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - all skylarks now migrate (#90)
 
+- `insectbiomass()` uses units
 
 ## [v0.5.1] - 13-06-2024
 
diff --git a/src/crop/farmplot.jl b/src/crop/farmplot.jl
index 13b7ad3..596da80 100644
--- a/src/crop/farmplot.jl
+++ b/src/crop/farmplot.jl
@@ -196,7 +196,7 @@ Return the height of the crop at this position, or nothing if there is no crop h
 (utility wrapper).
 """
 function cropheight(pos::Tuple{Int64,Int64}, model::SimulationModel)
-    ismissing(model.landscape[pos...].fieldid) ? nothing :
+    ismissing(model.landscape[pos...].fieldid) ? 0cm : #FIXME should not return 0
               model.farmplots[model.landscape[pos...].fieldid].height
 end
 
@@ -208,7 +208,7 @@ here (utility wrapper).
 """
 function cropcover(pos::Tuple{Int64,Int64}, model::SimulationModel)
     #FIXME LAItotal != ground cover?
-    ismissing(model.landscape[pos...].fieldid) ? nothing :
+    ismissing(model.landscape[pos...].fieldid) ? 0 : #FIXME should not return 0
               model.farmplots[model.landscape[pos...].fieldid].LAItotal
 end
 
diff --git a/src/nature/individuals.jl b/src/nature/individuals.jl
index f1f2488..5a056e7 100644
--- a/src/nature/individuals.jl
+++ b/src/nature/individuals.jl
@@ -81,10 +81,10 @@ Add the given location to the animal's territory.
 """
 function occupy!(animal::Animal, model::SimulationModel, position::Tuple{Int64,Int64})
     if isoccupied(model, speciesof(animal), position) #XXX should this be an error?
-        @warn "Position $position is already occupied by a $(speciesof(animal))."
+        @warn "Position $position is already occupied by a $(speciesof(animal))." animal
     end
-    push!(animal.territories, position)
-    push!(model.landscape[position...], animal.id)
+    push!(animal.territory, position)
+    push!(model.landscape[position...].territories, animal.id)
 end
 
 """
@@ -144,7 +144,10 @@ Let the animal move a given number of steps in the given direction ("north", "no
 "east", "southeast", "south", "southwest", "west", "northwest", "random").
 """
 function walk!(animal::Animal, model::SimulationModel, direction::String, distance::Length=-1m)
-    distance < 0m ? steps = 1 : steps = Int(floor(distance/@param(world.mapresolution)))
+    steps = 1
+    if distance > @param(world.mapresolution)
+        steps = Int(floor(distance/@param(world.mapresolution)))
+    end
     if direction == "north"
         shift = (0,-steps)
     elseif direction == "northeast"
@@ -162,7 +165,7 @@ function walk!(animal::Animal, model::SimulationModel, direction::String, distan
     elseif direction == "northwest"
         shift = (-steps,-steps)
     elseif direction == "random"
-        shift = Tuple(@rand([-steps:steps], 2))
+        shift = Tuple(@rand(-steps:steps, 2))
     else
         @error "Invalid direction in @walk: "*direction
     end
diff --git a/src/nature/macros.jl b/src/nature/macros.jl
index 1fa21f4..add747c 100644
--- a/src/nature/macros.jl
+++ b/src/nature/macros.jl
@@ -479,7 +479,6 @@ used nested within [`@phase`](@ref).
 macro walk(args...)
     #XXX add `ifempty` keyword?
     :(walk!($(esc(:self)), $(esc(:model)), $(map(esc, args)...)))
-    #FIXME MethodError: no method matching walk!(::Persefone.Skylark, ::AgricultureModel, ::Tuple{UnitRange{Int64}, UnitRange{Int64}}) -> due to map(esc())?
 end
 #TODO add own walking functions that respect habitat descriptors
 
diff --git a/src/nature/nature.jl b/src/nature/nature.jl
index d236343..317a19f 100644
--- a/src/nature/nature.jl
+++ b/src/nature/nature.jl
@@ -37,7 +37,7 @@ function speciesof(a::Union{Animal,Type})
     # strip out the module name if necessary (`Persefone.<species>`)
     (a isa Animal) && (a = typeof(a))
     spstrings = split(string(a), ".")
-    length(spstrings) == 1 ? spstrings[1] : spstrings[2]
+    length(spstrings) == 1 ? String(spstrings[1]) : String(spstrings[2])
 end
 
 """
diff --git a/src/nature/populations.jl b/src/nature/populations.jl
index e47bcd2..d48d759 100644
--- a/src/nature/populations.jl
+++ b/src/nature/populations.jl
@@ -145,8 +145,8 @@ end
 Test whether this location is part of the territory of an animal of the given species.
 """
 function isoccupied(model::SimulationModel, species::String, position::Tuple{Int64,Int64})
-    for terr in model.landscape[position...].territories
-        (speciesof(model.animals[terr]) == species) && return true
+    for id in model.landscape[position...].territories
+        (speciesof(model.animals[id]) == species) && return true
     end
     return false
 end
diff --git a/src/nature/species/skylark.jl b/src/nature/species/skylark.jl
index c564205..fd366ab 100644
--- a/src/nature/species/skylark.jl
+++ b/src/nature/species/skylark.jl
@@ -47,6 +47,7 @@ At the moment, this implementation is still in development.
     const migrationmortality::Float64 = 0.33 # chance of dying during the winter
 
     const minimumterritory = 5000m² # size of territory under ideal conditions
+    const mindistancetoedge = 60m # minimum distance of habitat to vertical structures
     const maxforageheight = 50cm # maximum preferred vegetation height for foraging
     const maxforagecover = 0.7 # maximum preferred vegetation cover for foraging
     const nestingheight = (15cm, 25cm) # min and max preferred vegetation height for nesting
@@ -101,16 +102,20 @@ end
 Males returning from migration move around to look for suitable habitats to establish a territory.
 """
 @phase Skylark territorysearch begin
-    #TODO
-    #TODO Standorttreue
-    # If we've found a territory, or the breeding season is over, move to the next phase
     if !isempty(self.territory)
-        @setphase(occupation)
+        @setphase(occupation) # If the male still has a territory, it is occupied again
     elseif month(model.date) > self.nestingend
-        @setphase(nonbreeding)
-    else
-        @walk("random", self.movementrange)
-        #FIXME MethodError: no method matching walk!(::Persefone.Skylark, ::AgricultureModel, ::Tuple{UnitRange{Int64}, UnitRange{Int64}})
+        @setphase(nonbreeding) # If the breeding season is over, give up looking
+    else # otherwise check whether this is a suitable location
+        newterritory = findterritory(self, model)
+        if isempty(newterritory)
+            @walk("random", self.movementrange)
+        else
+            for p in newterritory
+                @occupy(p)
+            end
+            @debug("$(animalid(self)) has found a territory.")
+        end
     end
 end
 
@@ -121,7 +126,10 @@ adjusting it to new conditions when and as necessary.
 @phase Skylark occupation begin
     #move to a random location in the territory
     @move(@rand(self.territory))
-    #TODO adjust territory as needed
+    if month(model.date) > self.nestingend #XXX check that all the young have left the nest
+        @setphase(nonbreeding)
+    end
+    #TODO adjust territory as needed (e.g. once a week, or when a brood is done?)
 end
 
 """
@@ -157,36 +165,35 @@ Females that have found a partner build a nest and lay eggs in a suitable locati
     if model.date < Date(year(model.date), self.nestingbegin...)
         # wait for nesting to begin, moving around in the territory
         @move(@rand(@animal(self.mate).territory))
-    elseif isempty(nest)
+    elseif month(model.date) > self.nestingend
+        @setphase(nonbreeding) # stop trying to nest if it's too late (should be rare)
+    elseif isempty(self.nest)
         # choose site, build nest & lay eggs
         for pos in @shuffle!(deepcopy(@animal(self.mate).territory))
-            #TODO is this condition correct? -> needs validation!
-            if (@landcover() == grass || @landcover() == soil ||
-                (@landcover() == agriculture &&
-                 (self.nestingheight[1] <= @cropheight() <= self.nestingheight[2]) &&
-                 (self.nestingcover[1] <= @cropcover() <= self.nestingcover[2])))
+            if allowsnesting(self, model, pos)
                 @move(pos)
                 self.nest = pos
                 self.clutch = @rand(self.eggsperclutch)
-                timer = self.nestbuildingtime + self.clutch # time to build + 1 day per egg laid
+                # time to build + 1 day per egg laid
+                self.timer = @rand(self.nestbuildingtime) + self.clutch
                 if month(model.date) == self.nestingbegin[1] 
                     # the first nest takes twice as long to build
                     #XXX this may affect the first two nests
-                    timer += self.nestbuildingtime
+                    self.timer += self.timer*2-self.clutch
                 end
                 break
             end
         end
-        isempty(nest) && @warn("$(animalid(self)) didn't find a nesting location.")
-    elseif timer == 0
+        isempty(self.nest) && @warn("$(animalid(self)) didn't find a nesting location.")
+    elseif self.timer == 0
         @debug("$(animalid(self)) has laid $(self.clutch) eggs.")
         @setphase(breeding)
     else
         self.timer -= 1
     end
     # tillage and harvest destroys the nest
-    @respond(tillage, nest = ())
-    @respond(harvesting, nest = ())
+    @respond(tillage, self.nest = ())
+    @respond(harvesting, self.nest = ())
 end
 
 """
@@ -196,7 +203,7 @@ chicks are independent or in case of brood loss.
 @phase Skylark breeding begin
     #TODO wait for eggs to hatch & chicks to mature, checking for mortality
     # restart breeding cycle if there is time
-    if clutch == 0 && month(model.date) <= self.nestingend
+    if self.clutch == 0 && month(model.date) <= self.nestingend
         @setphase(nesting) #TODO breeding delay?
     elseif month(model.date) > self.nestingend
         @setphase(nonbreeding)
@@ -212,7 +219,7 @@ end
 Select the dates on which this skylark will leave for / return from its migration,
 based on observed migration patterns.
 """
-function migrationdates(skylark::Animal, model::SimulationModel)
+function migrationdates(skylark::Skylark, model::SimulationModel)
     #TODO this ought to be temperature-dependent and dynamic
     minleave = skylark.sex == female ? (September, 15) : (October, 1)
     minarrive = skylark.sex == male ? (February, 15) : (March, 1)
@@ -224,29 +231,95 @@ function migrationdates(skylark::Animal, model::SimulationModel)
 end
 
 """
-    foragequality(maxcover, maxheight, position, model)
+    findterritory(skylark, model)
+
+Check whether the habitat surrounding the skylark is suitable for establishing a territory.
+If it is, return the list of coordinates that make up the new territory, else return an empty list.
+"""
+function findterritory(skylark::Skylark, model::SimulationModel)
+    effectivesize::Area = 0m² # the usable size of the territory, weighted by habitat quality
+    territory::Vector{Tuple{Int64,Int64}} = []
+    msize = size(model.landscape)
+    radius = 0
+    constrained = false
+    # Inspect the landscape in concentric circles around the individual until enough pixels have
+    # been found to provide a territory of sufficient size and quality. If there are no suitable
+    # pixels in one circle, break off the search (territories must be contiguous).
+    while !constrained && effectivesize < skylark.minimumterritory
+        constrained = true
+        if radius == 0
+            coords = [skylark.pos]
+        else # list all coordinates in the next circle...
+            coords = []
+            xrange = (skylark.pos[1]-radius, skylark.pos[1]+radius)
+            yrange = (skylark.pos[2]-radius, skylark.pos[2]+radius)
+            for x in xrange[1]:xrange[2]
+                push!(coords, (x, yrange[1]))
+                push!(coords, (x, yrange[2]))
+            end
+            for y in (yrange[1]+1):(yrange[2]-1) #avoid duplicating the corners
+                push!(coords, (y, xrange[1]))
+                push!(coords, (y, xrange[2]))
+            end
+        end
+        #FIXME some duplicates remain?
+        for c in coords # ...then inspect them
+            (c[1] <= 0 || c[2] <= 0 || c[1] > msize[1] || c[2] > msize[2]) && continue
+            (isoccupied(model, "Skylark", c)) && continue
+            quality = foragequality(skylark, model, c)
+            #XXX check for nesting habitats?
+            if quality > 0
+                constrained = false
+                push!(territory, c)
+                effectivesize += @param(world.mapresolution)^2*quality
+            end
+        end
+        radius +=1
+    end
+    constrained ? [] : territory
+end
+
+"""
+    foragequality(skylark, model, pos)
 
 Calculate the relative quality of the habitat at this position for foraging.
 (Approximated from Püttmanns et al., 2021; Jeromin, 2002; Jenny, 1990b.)
 """
-function foragequality(maxcover::Float64, maxheight::Float64, position::Tuple{Int64,Int64},
-                       model::SimulationModel)
+function foragequality(skylark::Skylark, model::SimulationModel, pos::Tuple{Int64,Int64})
     #TODO this is a key function that needs to be validated thoroughly
-    px = model.landscape[position...]
-    !(px.landcover in (grass, soil, agriculture)) && return 0.0
+    if !(@landcover() in (grass, soil, agriculture)) ||
+        (@distanceto(forest) < skylark.mindistancetoedge) ||
+        (@distanceto(builtup) < skylark.mindistancetoedge)
+        return 0.0
+    end
     quality = 1.0
-    if !is.missing(px.fieldid)
-        f = model.farmplots[px.fieldid]
-        groundcoverfactor = x -> bounds(-1/maxcover + 1.0, max=1.0)
-        plantheightfactor = x -> bounds(-1/maxheight + 1.0, max=1.0)
-        #FIXME need percentage cover, not LAI
-        quality = bounds(groundcoverfactor(f.LAItotal) + plantheightfactor(f.height), max=1.0)
-    else
-        @warn "No field assigned to location $position" px
+    f = farmplot(pos, model)
+    # Assume that grass and soil have a habitat quality of 1.0. For fields, calculate quality
+    # as the sum of cover quality and height quality (each modelled using a linear regression).
+    if !isnothing(f)
+        groundcoverfactor = x -> bounds((-1/skylark.maxforagecover)*x + 1.0, max=1.0)
+        plantheightfactor = x -> bounds((-1/skylark.maxforageheight)*(x |> cm) + 1.0, max=1.0)
+        #FIXME need percentage cover in FarmPlot, not LAI
+        #FIXME height is currently dimensionless in FarmPlot, hence the conversion to cm
+        quality = bounds(groundcoverfactor(f.LAItotal) + plantheightfactor(f.height*1cm), max=1.0)
     end
     return quality
 end
 
+"""
+    allowsnesting(skylark, model, pos)
+
+Check whether the given position is suitable for nesting.
+"""
+function allowsnesting(skylark::Skylark, model::SimulationModel, pos::Tuple{Int64,Int64})
+    #TODO is this condition correct? -> needs validation!
+    (@landcover() == grass ||
+     (@landcover() == agriculture &&
+      (skylark.nestingheight[1] <= @cropheight() <= skylark.nestingheight[2]) &&
+      (skylark.nestingcover[1] <= @cropcover() <= skylark.nestingcover[2]))) #&&
+      #(@distanceto(forest) < skylark.mindistancetoedge) &&
+      #(@distanceto(builtup) < skylark.mindistancetoedge)
+end
 
 ## INITIALISATION
 
@@ -274,9 +347,9 @@ end
 @populate Skylark begin
     # initialise on open land, at least 60m from vertical structures
     habitat = @habitat((@landcover() == grass || @landcover() == agriculture) &&
-                       @distanceto(forest) >= 60m && @distanceto(builtup) >=60m)
+                       @distanceto(forest) >= 60m && @distanceto(builtup) >=60m) #XXX magic numbers
     initphase = nonbreeding
     birthphase = nonbreeding
-    indarea = 3ha
+    indarea = 3ha #XXX magic number
     pairs = true
 end
-- 
GitLab