Incorporate density-dependent growth

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: FLBEIA
 Title: Bio-Economic Impact Assessment of Management Strategies using FLR
-Version: 1.16.1.11
-Date: 2023-03-07
+Version: 1.16.1.13
+Date: 22024-02-16
 Authors@R: c(person("FLBEIA", "Team", email = "[email protected]", role = c("aut","cre")))
 Description: A simulation toolbox that describes a fishery system under
   a Management Strategy Estrategy approach. The objective of the model is
@@ -121,6 +121,7 @@ Collate:
     "OM_1a_Population_Growth_Functions.R"
     "OM_1a1_Stock_Recruitment_functions.R"
     "OM_1a1_Stock_Recruitment_functions_segregmix.R"
+    "OM_1a2_DensityDependent_weight_functions.R"
     "OM_2_fleet.om.R"
     "OM_2a_Effort_Dynamics.R"
     "OM_2a_Effort_Dynamics_Auxiliary.R"

R/Conditioning_create.ecoData.arrays.R

@@ -137,7 +137,7 @@ create.ecoData <- function(file, fltObj, hist.yrs, mean.yrs, sim.yrs){
     covars[["w2"]][fl,ac(w2$year)]                   <- w2[,'fleet']
     covars[["NewVesselCost"]][fl,ac(nev$year)]       <- nev[,'fleet']
     covars[["EmploymentPerVessel"]][fl,ac(emp$year)] <- emp[,'fleet']
-    covars[["Depreciation"]][fl,ac(emp$year)]        <- dep[,'fleet']
+    covars[["Depreciation"]][fl,ac(dep$year)]        <- dep[,'fleet']
   }
 
     #---------------------------------------------------------------------------

R/MAIN_Function.R

@@ -369,6 +369,7 @@ FLBEIA <- function(biols, SRs = NULL, BDs = NULL, fleets, covars = NULL, indices
         biols <- res$biols
         SRs   <- res$SRs
         BDs   <- res$BDs
+        covars <- res$covars
 
         cat('------------ FLEETS OM ------------\n')
         # - Fleets OM.

R/MP_2_Assessment_spict2flbeia.R

@@ -75,9 +75,9 @@ spict2flbeia <- function(stock, indices, control=NULL,covars=covars){
 
   }
 
-  stock@stock[,years] <- results[[st]][["spict_B"]] #explotable biomass
+  stock@stock[,years] <- results[[st]][["spict_B"]] #exploitable biomass
   stock@stock.n[,years] <- stock@stock/stock@stock.wt
-  stock@harvest[,years] <- results[[st]][["spict_F"]] #explotable biomass
+  stock@harvest[,years] <- results[[st]][["spict_F"]] #fishing mortality
   covars[[st]]$spict_Bmsy[,tail(years,n=1)] <- results[[st]][["spict_Bmsy"]][,tail(years,n=1)]
   covars[[st]]$spict_Fmsy[,tail(years,n=1)] <- results[[st]][["spict_Fmsy"]][,tail(years,n=1)]
   # Returns what?

R/OM_1_biol.om.R

@@ -27,11 +27,13 @@ biols.om <- function(biols, fleets, SRs, BDs, covars, biols.ctrl, year, season){
 
         biols[[st]] <- res$biol
 
+        if(!is.null(res$covars)) covars <- res$covars
+
         if(!is.null(SRs[[st]]))  SRs[[st]] <- res$SR
         if(!is.null(BDs[[st]]))  BDs[[st]] <- res$BD
     }
 
-    return(list(biols = biols, SRs = SRs, BDs = BDs))
+    return(list(biols = biols, SRs = SRs, BDs = BDs, covars = covars))
 
 }
 

R/OM_1a2_DensityDependent_weight_functions.R

@@ -0,0 +1,63 @@
+#-------------------------------------------------------------------------------
+#                   DENSITY-DEPENDENT WEIGHT FUNCTIONS
+#
+# 14/02/2024
+#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# REMARK: '...' in the arguments of the functions are necessary in order to be
+#   generalist inside 'ASPG_DDW' function ('eval(call(...)').
+#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+# Input variables from ASPG_DWW call:
+# biol = biol, stknm = stknm, year = year, season = season, ctrl = ddw.ctrl, covars = covars
+
+
+#-------------------------------------------------------------------------------
+# ddwAgeCa(biol, SR, fleets, biol.control)
+# - OUTPUT: list(wt = wage, wt.chg = wt.chg) - vector with estimated weight at age values and relative change
+#-------------------------------------------------------------------------------
+
+# Weights-at-age based on SSB and a covariate
+
+ddwAgeCa <- function(biol, stknm, year, season, ctrl, covars, ...) {
+
+  pars  <- ctrl$params[,,year,season,] # array[npar,nage,nyr,ns,nit]
+  covnm <- ctrl$covnm                  # name of the covariate
+
+  logssb  <- log(ssb(biol)[,year,,season,drop=TRUE])
+  wage    <- exp(pars["a",] + pars["b",] * logssb + pars["c",] * covars[[covnm]][stknm,year,,season,drop=TRUE])
+
+  # wt change
+  wt.ref <- biol@wt[,year,,season,]
+  wt.chg <- wage/wt.ref
+
+  return(list(wt = wage, wt.chg = wt.chg))
+
+}
+
+
+#-------------------------------------------------------------------------------
+# ddwAgeLFD(biol, SR, fleets, biol.control)
+# - OUTPUT: list(wt = wage, wt.chg = wt.chg) - vector with estimated weight at age values and relative change
+#-------------------------------------------------------------------------------
+
+# Weights-at-age based on SSB (linear model for estimating LW b parameter) and A LFD
+
+ddwAgeLFD <- function(biol, stknm, year, season, ctrl, covars, ...) {
+
+  lfd       <- ctrl[['LFD']]
+  a         <- ctrl[['a.lw']]
+  lbins     <- as.numeric(colNames(lfd))
+
+  B <- quantSums((biol@wt*biol@n)[,year-1])[drop=T] #! DG needs to consider season dimension
+
+  condF <- predict(LW_lm, data.frame(biomass = B))  #! DG requires: biols.ctrl[[stknm]][['ddw.ctrl']][['LW_lm']]
+
+  wy <- a*(lbins)^condF
+
+  wt. <- rowSums(sweep(lfd, 2, wy, "*"))  #! DG needs to consider also season dimension
+  wt  <- biol@wt[,year,,season,]
+
+  return(list(wt = wt., wt.chg = wt./wt))
+
+}
+

R/OM_1a_Population_Growth_Functions.R

@@ -32,7 +32,7 @@ fixedPopulation <- function(biols, SRs, fleets, year, season, stknm, ...)  retur
 
 #-------------------------------------------------------------------------------
 # ASPG(biol, SR, fleets, biol.control)
-# - OUTPUT: list(biol = biol, SR = SR) - Upadated FLBiol and FLSRsim objects.
+# - OUTPUT: list(biol = biol, SR = SR) - Updated FLBiol and FLSRsim objects.
 #-------------------------------------------------------------------------------
 
 ASPG <- function(biols, SRs, fleets, year, season, stknm, ...){
@@ -349,30 +349,38 @@ ASPG_Baranov <- function(biols, SRs, fleets, year, season, stknm, ...){
 } 
 
 
+#-------------------------------------------------------------------------------
+# ASPG_DDW(biol, SR, fleets, biol.control)
+# - OUTPUT: list(biol = biol, SR = SR, covars = covars) - Updated FLBiol and FLSRsim objects.
+#-------------------------------------------------------------------------------
+
 # Age structured population growth with densodependence
-#### ASPG_DDW ----
-ASPG_DDW <- function(biols, SRs, fleets, year, season, stknm, biols.ctrl,...){
-
-  # update wt of stknm
-  ddw.model <- biols.ctrl[[stkm]][['growth.model']]
-  lfd       <- biols.ctrl[[stkm]][['LFD']]
-  a         <- biols.ctrl[[stkm]][['a.lw']]
-  lbins     <- as.numeric(colNames(lfd))
+
+ASPG_DDW <- function(biols, SRs, fleets, stknm, year, season, ctrl, covars, ...){
 
-  B <- quantSums((biols[[stkm]]@wt*biols[[stkm]]@n)[,year-1])[drop=T]
+  biol <- biols[[stknm]]
 
-  condF <- predict(LW_lm, data.frame(biomass = B))
+  # check if DDW covars available
+  if (!"DDW" %in% names(covars) | !stknm %in% names(covars[["DDW"]]))
+    covars[["DDW"]][[stknm]] <- biols[[stknm]]@n * 0 + 1
 
-  wy <- a*(lbins)^condF
+  # update wt of stknm
+  ddw.model <- ctrl[[stknm]][['ddw.model']]
+  ddw.ctrl  <- ctrl[[stknm]][['ddw.ctrl']]
 
-  wt <- rowSums(sweep(lfd, 2, wy, "*"))
+  wts <- eval(call(ddw.model, biol = biol, stknm = stknm, year = year, season = season, 
+                  ctrl = ddw.ctrl, covars = covars))
 
-  biols[[stkm]]@wt[,year] <- wt
+  # save
+  biols[[stknm]]@wt[,year,,season,]        <- wts$wt
+  covars[["DDW"]][[stknm]][,year,,season,] <- wts$wt.chg
 
   # use normal ASPG to project the population
 
   res <- ASPG(biols, SRs, fleets, year, season, stknm, biols.ctrl,...)
 
+  res$covars <- covars
+
   return(res)
 }
 

R/OM_2_fleet.om.R

@@ -24,6 +24,20 @@ fleets.om <- function(fleets, biols, BDs, covars, advice, biols.ctrl, fleets.ctr
     fleets         <- unclass(fleets) # unclass fleets to speed up the algorithm
     fleets.ctrl.aux <- fleets.ctrl
 
+
+    # update weights if DDW function applied in any of the stocks
+    if(!(is.null(covars$DDW))) {
+      for(st in names(covars$DDW)){
+        for(fl in  names(which(sapply(lapply(fleets, catchNames), function(x) (st %in% x))))) {
+          for(mt in names(which(sapply(lapply(fleets[[fl]]@metiers, catchNames), function(x) (st %in% x))))) {
+            mult <- covars[['DDW']][[st]][,year,,season,]
+            fleets[[fl]]@metiers[[mt]]@catches[[st]]@landings.wt[,year,,season,] <- mult*fleets[[fl]]@metiers[[mt]]@catches[[st]]@landings.wt[,year,,season,]
+            fleets[[fl]]@metiers[[mt]]@catches[[st]]@discards.wt[,year,,season,] <- mult*fleets[[fl]]@metiers[[mt]]@catches[[st]]@discards.wt[,year,,season,]
+          }
+        }
+      }
+    }
+
     # 1. Calculate effort. 
     print('~~~~~~~~~~ EFFORT ~~~~~~~~')
 
@@ -44,7 +58,7 @@ fleets.om <- function(fleets, biols, BDs, covars, advice, biols.ctrl, fleets.ctr
                     year = year, season = season, biols.ctrl=biols.ctrl, fleets.ctrl = fleets.ctrl, covars = covars, 
                     assess.ctrl=assess.ctrl, advice.ctrl = advice.ctrl)) 
 
-        fleets[[fl]]         <- res$fleets[[fl]]
+        fleets[[fl]]          <- res$fleets[[fl]]
         fleets.ctrl.aux[[fl]] <- res$fleets.ctrl[[fl]]
         remove(res)
     }

README.md

@@ -1,6 +1,6 @@
 # FLBEIA
-- Version: 1.16.1.11
-- Date: 2023-03-07
+- Version: 1.16.1.13
+- Date: 2024-02-16
 - Author: Dorleta GARCIA <[email protected]>; FLBEIA Team <[email protected]>
 - Maintainer: Dorleta GARCIA, AZTI & FLBEIA Team
 - Repository: <https://github.com/flr/FLBEIA/>