Improve density-dependent growth function

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: FLBEIA
 Title: Bio-Economic Impact Assessment of Management Strategies using FLR
-Version: 1.16.1.13
-Date: 22024-02-16
+Version: 1.16.1.14
+Date: 22024-02-20
 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

NEWS

@@ -1,3 +1,11 @@
+#~~~~~~~~~~~~~~~~
+# 2024/02/20
+#~~~~~~~~~~~~~~~~
+FLBEIA 1.16.1.14
+- Include density-dependent growth function
+- Specific funcion for density-dependent weights for a stock with 4 age classes and linear relationship for weights estimation and abundance (dependent on SSB and SST).
+- Function using length-frequency distribution and previous year abundance is in progress.
+
 #~~~~~~~~~~~~~~~~
 # 2014/10/01
 #~~~~~~~~~~~~~~~~

R/OM_1a2_DensityDependent_weight_functions.R

@@ -17,17 +17,61 @@
 #-------------------------------------------------------------------------------
 
 # Weights-at-age based on SSB and a covariate
+# (specific for Bay of Biscay anchovy with 3 age classes - 0 to 3+)
 
 ddwAgeCa <- function(biol, stknm, year, season, ctrl, covars, ...) {
 
+  require(nleqslv)
+
   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.ref <- biol@wt[,year,,season,drop=TRUE]
+  ssna   <- (n(biol) * mat(biol) * exp(-spwn(biol) * m(biol)))[,year,,season,drop=TRUE] # exclude age0 (inmature)
+  sst    <- covars[[covnm]][stknm,year,,season,drop=TRUE]
+
+  # ssb.ref <- quantSums(ssna * wt.ref[-1])
+
+  wfun <- function(x) {
+
+    # x[i] is log(wage_i)
+    # x[i] = a_i + b_i * log(ssb) + c_i * sst
+    #      = a_i + b_i * log(sum_j[ssna_j * exp(x[j])]) + c_i * sst 
+    #         with ssna_j = @n_j * @mat_j * exp(-@m_j * @spwn_j)
+
+    y <- numeric(3)
+
+    logssb <- log(ssna[1] * exp(x[1]) + ssna[2] * exp(x[2]) +
+                    ssna[3] * exp(x[3]) + ssna[4] * exp(x[4]))
+
+    y[1] <- x[1] - pars["b",1] * logssb - pars["a",1] - pars["c",1] * sst
+    y[2] <- x[2] - pars["b",2] * logssb - pars["a",2] - pars["c",2] * sst
+    y[3] <- x[3] - pars["b",3] * logssb - pars["a",3] - pars["c",3] * sst
+    y[4] <- x[4] - pars["b",4] * logssb - pars["a",4] - pars["c",4] * sst
+
+    return(y)
+  }
+
+  lwest <- nleqslv(log(wt.ref), fn = wfun, control = list(btol=1e-08, delta="newton"))
+
+  if (lwest$termcd != 1)
+    print(paste0(stknm," weights calculation message: ",lwest3$message))
+
+  if (any(wfun(lwest$x)>1e-08))
+    stop(paste("Issues with mean weight calculation for", stknm))
+
+  wage <- exp(lwest$x)
+
+  # check
+  logssb <- sum(ssna * wage)
+  di <- numeric(4)
+  for (i in 1:4) {
+    di[i] <- wage[i] - exp(pars["a",i] + pars["b",i] * logssb + pars["c",i] * sst)
+  }
+  if (any(round(di-wage,8) != 0))
+    stop(paste("Issues with density-dependent weights-at-age for",stknm))
 
   # wt change
-  wt.ref <- biol@wt[,year,,season,]
   wt.chg <- wage/wt.ref
 
   return(list(wt = wage, wt.chg = wt.chg))

R/OM_1a_Population_Growth_Functions.R

@@ -358,16 +358,19 @@ ASPG_Baranov <- function(biols, SRs, fleets, year, season, stknm, ...){
 
 ASPG_DDW <- function(biols, SRs, fleets, stknm, year, season, ctrl, covars, ...){
 
-  biol <- biols[[stknm]]
+  cat('-----------------ASPG_DDW-----------\n')
 
   # check if DDW covars available
   if (!"DDW" %in% names(covars) | !stknm %in% names(covars[["DDW"]]))
     covars[["DDW"]][[stknm]] <- biols[[stknm]]@n * 0 + 1
 
-  # update wt of stknm
+  # estimate the DD weights
   ddw.model <- ctrl[[stknm]][['ddw.model']]
   ddw.ctrl  <- ctrl[[stknm]][['ddw.ctrl']]
 
+  if (ddw.model == "ddwAgeCa") # for in-year ssb calculation (nage required)
+    biol <- ASPG(biols, SRs, fleets, year, season, stknm, biols.ctrl,...)$biol
+
   wts <- eval(call(ddw.model, biol = biol, stknm = stknm, year = year, season = season, 
                   ctrl = ddw.ctrl, covars = covars))
 

README.md

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