From 3936481dc2b2ad1150cba2aed2a2dc7e8afe5863 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sonia=20S=C3=A1nchez?= Date: Wed, 17 Jul 2024 09:15:08 +0200 Subject: [PATCH] Adapt ddwAgeCa to several iterations --- DESCRIPTION | 4 +- NEWS | 15 ++++ R/OM_1a2_DensityDependent_weight_functions.R | 85 +++++++++++++------- README.md | 4 +- 4 files changed, 77 insertions(+), 31 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index df1266d..9673f9f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: FLBEIA Title: Bio-Economic Impact Assessment of Management Strategies using FLR -Version: 1.16.1.16 -Date: 2024-06-28 +Version: 1.16.1.19 +Date: 2024-07-16 Authors@R: c(person("FLBEIA", "Team", email = "flbeia@azti.es", 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 diff --git a/NEWS b/NEWS index 2622bbb..901d0b8 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,18 @@ +#~~~~~~~~~~~~~~~~ +# 2024/07/12 +#~~~~~~~~~~~~~~~~ +FLBEIA 1.16.1.19 +- XXXXX + + + +#~~~~~~~~~~~~~~~~ +# 2024/07/09 +#~~~~~~~~~~~~~~~~ +FLBEIA 1.16.1.17 +- New price functions available (elasticPriceAge) and previous ones (elasticPrice) include alternative options now. + + #~~~~~~~~~~~~~~~~ # 2024/02/20 #~~~~~~~~~~~~~~~~ diff --git a/R/OM_1a2_DensityDependent_weight_functions.R b/R/OM_1a2_DensityDependent_weight_functions.R index 9481c33..92a03c8 100644 --- a/R/OM_1a2_DensityDependent_weight_functions.R +++ b/R/OM_1a2_DensityDependent_weight_functions.R @@ -23,16 +23,21 @@ ddwAgeCa <- function(biol, stknm, year, season, ctrl, covars, ...) { require(nleqslv) - pars <- ctrl$params[,,year,season,] # array[npar,nage,nyr,ns,nit] + na <- dim(biol@wt)[1] + it <- dim(biol@wt)[6] + + pars <- ctrl$params[,,year,season,] # array[npar,nage,1,1,nit] covnm <- ctrl$covnm # name of the covariate - 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] + wt.ref <- biol@wt[,year,,season,drop=TRUE] # [na,nit] or [na] + ssna <- (n(biol) * mat(biol) * exp(-spwn(biol) * m(biol)))[,year,,season,drop=TRUE] # [na,nit] or [na] # to exclude age0 (immature) + sst <- covars[[covnm]][stknm,year,,season,drop=TRUE] # [nit] + + wage <- wt.ref * NA # ssb.ref <- quantSums(ssna * wt.ref[-1]) - wfun <- function(x) { + wfun <- function(x, param, ssna) { # x[i] is log(wage_i) # x[i] = a_i + b_i * log(ssb) + c_i * sst @@ -44,33 +49,59 @@ ddwAgeCa <- function(biol, stknm, year, season, ctrl, covars, ...) { 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 + y[1] <- x[1] - param["b",1] * logssb - param["a",1] - param["c",1] * sst + y[2] <- x[2] - param["b",2] * logssb - param["a",2] - param["c",2] * sst + y[3] <- x[3] - param["b",3] * logssb - param["a",3] - param["c",3] * sst + y[4] <- x[4] - param["b",4] * logssb - param["a",4] - param["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 (it == 1) { + + lwest <- nleqslv(log(wt.ref), fn = wfun, param=pars, ssna=ssna, control = list(btol=1e-08, delta="newton")) + + if (lwest$termcd != 1) + print(paste0(stknm," weights calculation message: ",lwest3$message)) + + if (any(wfun(lwest$x, param=pars, ssna=ssna)>1e-08)) + stop(paste("Issues with mean weight calculation for", stknm)) + + wage <- exp(lwest$x) + + # # check + # logssb <- sum(ssna * wage) + # di <- numeric(na) + # for (a in 1:na) + # di[a] <- wage[a] - exp(pars["a",a] + pars["b",a] * logssb + pars["c",a] * sst) + # if (any(round(di-wage,4) != 0)) + # stop(paste("Issues with density-dependent weights-at-age for",stknm)) + + } else { + + for (i in 1:it) { + + lwest <- nleqslv(log(wt.ref[,i]), fn = wfun, param=pars[,,i], ssna=ssna[,i], control = list(btol=1e-08, delta="newton")) + + if (lwest$termcd != 1) + print(paste0(stknm," weights calculation message: ",lwest3$message)) + + if (any(wfun(lwest$x, param=pars[,,i], ssna=ssna[,i])>1e-08)) + stop(paste("Issues with mean weight calculation for", stknm)) + + wage[,i] <- exp(lwest$x) + } + + # # check + # logssb <- apply(ssna * wage, 2, sum) + # di <- array(NA, dim=c(na,it)) + # for (a in 1:na) + # di[a,] <- wage[a,] - exp(pars["a",a,] + pars["b",a,] * logssb + pars["c",a,] * sst) + # if (any(round(di-wage,4) != 0)) + # stop(paste("Issues with density-dependent weights-at-age for",stknm)) + } - if (any(round(di-wage,8) != 0)) - stop(paste("Issues with density-dependent weights-at-age for",stknm)) - + # wt change wt.chg <- wage/wt.ref diff --git a/README.md b/README.md index 5e8b241..6c29a7e 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # FLBEIA -- Version: 1.16.1.18 -- Date: 2024-07-09 +- Version: 1.16.1.19 +- Date: 2024-07-16 - Author: Dorleta GARCIA ; FLBEIA Team - Maintainer: Dorleta GARCIA, AZTI & FLBEIA Team - Repository: