#' Reads bootstrap estimates from a stock assessment model (including VPA) into an operating model. \code{Assess2OM}
#' is identical to \code{VPA2OM}.
#'
#'
#' @description A function that uses a set of bootstrap estimates of numbers-at-age, fishing mortality rate-at-age, M-at-age,
#' weight-at-age, length-at-age and Maturity-at-age to define a fully described MSEtool operating model. The user still
#' needs to parameterize most of the observation and implementation portions of the operating model.
#' @param Name Character string. The name of the operating model.
#' @param proyears Positive integer. The number of projection years for MSE.
#' @param interval Positive integer. The interval at which management procedures will update the management advice in \link[MSEtool]{runMSE}, e.g., 1 = annual updates.
#' @param CurrentYr Positive integer. The current year (final year of fitting to data)
#' @param h The steepness of the stock-recruitment curve (greater than 0.2 and less than 1, assumed to be close to 1 to match VPA assumption).
#' Either a single numeric or a length nsim vector.
#' @param Obs The observation model (class Obs). This function only updates the catch and index observation error.
#' @param Imp The implementation model (class Imp). This function does not update implementation parameters.
#' @param naa Numeric array `[sim, ages, year]`. Numbers-at-age `[first age is age zero]`.
#' @param faa Numeric array `[sim, ages, year]`. Fishing mortality rate-at-age `[first age is age zero]`.
#' @param waa Numeric array `[sim, ages, year]`. Weight-at-age `[first age is age zero]`.
#' @param Mataa Numeric array `[sim, ages, year]`. Maturity (spawning fraction)-at-age `[first age is age zero]`.
#' @param Maa Numeric array `[sim, ages, year]`. Natural mortality rate-at-age `[first age is age zero]`.
#' @param laa Numeric array `[sim, ages, year]`. Length-at-age `[first age is age zero]`.
#' @param nyr_par_mu Positive integer. The number of recent years that natural mortality, age vulnerability, weight, length and maturity parameters are averaged over for defining future projection conditions.
#' @param LowerTri Integer. The number of recent years for which model estimates of recruitment are ignored (not reliably estimated by the assessment)
#' @param recind Positive integer. The first age class that fish 'recruit to the fishery'. The default is 0 - ie the first position in the age dimension of naa is age zero
#' @param plusgroup Logical. Does the assessment assume that the oldest age class is a plusgroup?
#' @param altinit Integer. Various assumptions for how to set up the initial numbers. 0: standard, 1: no plus group, 2: temporary fix for MSEtool plus group initialization
#' @param fixq1 Logical. Should q be fixed (ie assume the F-at-age array faa is accurate?
#' @param report Logical, if TRUE, a diagnostic will be reported showing the matching of the OM reconstructed numbers at age vs the assessment.
#' @param silent Whether to silence messages to the console.
#' @param ... Additional arguments (for all, either a numeric or a length nsim vector):
#' \itemize{
#' \item \code{fecaa} Fecundity at age. Default fecundity is the product of maturity and weight at age.
#' \item \code{SRrel} Stock-recruit relationship. (\code{1} for Beverton-Holt (default), \code{2} for Ricker)
#' \item \code{R0} unfished recruitment
#' \item \code{phi0} unfished spawners per recruit associated with R0 and h. With time-varying parameters, openMSE uses the mean phi0
#' in the first \code{ageM} (age of 50 percent maturity) years for the stock-recruit relationship. \code{Assess2OM} will re-calculate R0 and h
#' in the operating model such that the stock-recruit \code{alpha} and \code{beta} parameters match values implied in the input.
#' \item \code{Perr} recruitment standard deviation (lognormal distribution) for sampling future recruitment
#' \item \code{AC} autocorrelation in future recruitment deviates.
#' \item \code{spawn_time_frac} The fraction of a year when spawning takes place (e.g., 0.5 is the midpoint of the year)
#' }
#' @details Use a seed for the random number generator to sample future recruitment.
#' @return An object of class \linkS4class{OM}.
#' @author T. Carruthers
#' @aliases VPA2OM
#' @export
#' @seealso \link{SS2OM} \link{iSCAM2OM} \link{WHAM2OM} \link{ASAP2OM}
Assess2OM <- function(Name="A fishery made by VPA2OM",
proyears=50, interval=2, CurrentYr=as.numeric(format(Sys.Date(), "%Y")),
h=0.999,
Obs = MSEtool::Imprecise_Unbiased, Imp=MSEtool::Perfect_Imp,
naa, faa, waa, Mataa, Maa, laa,
nyr_par_mu = 3, LowerTri=1,
recind=0, plusgroup=TRUE, altinit=0, fixq1=TRUE,
report=FALSE, silent=FALSE, ...) {
simup<-function(param, OM, do_slot = TRUE){ # Generic function for converting VPA outputs to an OM
paramnam<-deparse(substitute(param))
if(length(param)==1){
OM@cpars[[paramnam]]<-rep(param,nsim)
}else if(length(param)<nsim){
stop(paste('parameter vector', paramnam, 'was length',length(param),'which is neither length 1 or nsim long'))
}else{
OM@cpars[[paramnam]]<-param
}
if(do_slot) slot(OM,paramnam)<-rep(param[1],2)
OM
}
dots <- list(...) #R0, phi0
cond <- vapply(list(faa, waa, Mataa, Maa, laa), function(x) all(dim(naa) == dim(x)), logical(1))
if(!length(cond) || any(is.na(cond)) || any(!cond)) {
stop('One or more of the following arrays do not have the same shape: naa, faa, waa, Mataa, Maa, Laa')
}
if(!is.null(dots$fecaa)) {
fecaa <- dots$fecaa
if(!all(dim(naa) == dim(fecaa))) stop("Dimension of fecaa not equal to that for naa")
} else {
fecaa <- waa * Mataa
}
nyears<-dim(naa)[3]
nsim<-dim(naa)[1]
if(recind==1) { # create a dummy age 0 dimension to the various arrays
ageind<-1:dim(naa)[2]
dims<-c(nsim,1,nyears)
zeros<-array(0,dims)
# N0 back imputed
rec<-naa[,1,2:nyears]#*exp(Maa[,1,2:nyears])
muR0<-apply(rec,1,mean) # mean R0 is assumed for most recent N0 - will be filled anyway using LowerTri argument
N0<-array(cbind(rec,muR0),dims)
naa<-abind(N0,naa,along=2)
# F, weight, length and maturity assumed to be zero
faa<-abind(zeros,faa,along=2)
waa<-abind(zeros,waa,along=2)
Mataa<-abind(zeros,Mataa,along=2)
laa<-abind(zeros,laa,along=2)
fecaa<-abind(zeros,fecaa,along=2)
# M copied from first year to age zero
Maa<-abind(array(tiny + .Machine$double.eps,dim(Maa[,1,])),Maa,along=2)
message("Age zero positions for arrays were created with the following assumptions: N(0) = N(1) * exp(M(1)), N0 in most recent year is mean(R0), F(0) = weight(0) = maturity(0) = length(0) = 0, M(0) = M(1)")
}
# Dimensions
n_age <- dim(naa)[2]
maxage<-dim(naa)[2] - 1
# Stock-recruit relationship
if(!is.null(dots$SRrel)) {
SRrel <- dots$SRrel
} else {
SRrel <- 1
}
# Mine values to spec-out the OM
if(!is.null(dots$R0)) {
if(length(dots$R0) == 1) {
R0 <- rep(dots$R0, nsim)
} else {
R0 <- dots$R0
}
} else {
R0 <- naa[,1,1] # R0 is mean numbers at age 0
}
if(length(h) == 1) h <- rep(h, nsim)
OM<-new('OM')
if(!is.null(dots$spawn_time_frac)) {
if(lengths(dots$spawn_time_frac) == 1) {
spawn_time_frac <- rep(dots$spawn_time_frac, nsim)
} else {
spawn_time_frac <- dots$spawn_time_frac
}
} else {
spawn_time_frac <- rep(0, nsim)
}
SBsurv <- lapply(1:nsim, calc_survival,
StockPars = list(M_ageArray = Maa, n_age = n_age, spawn_time_frac = spawn_time_frac),
plusgroup = plusgroup, inc_spawn_time = TRUE) %>%
abind(along = 3) %>%
aperm(c(3,1,2))
SSBpR <- apply(SBsurv * fecaa, c(1, 3), sum)
SSBpR_out <- vapply(1:nsim, function(x) {
ageM <- min(LinInterp(Mataa[x,, 1], y = 1:n_age, 0.5), maxage - 1)
SSBpR[x, 1:(ceiling(ageM) + 1)] %>% mean()
}, numeric(1))
# Recalculate R0 and h that uses the replacement line corresponding to 1/SSBpR_out (used by openMSE)
if(!is.null(dots$phi0)) {
new_SR <- local({
if(length(dots$phi0) == 1) {
phi0 <- rep(dots$phi0, nsim)
} else {
phi0 <- dots$phi0
}
if(SRrel == 1) {
Arec <- 4*h/(1-h)/phi0
Brec <- (5*h-1)/(1-h)/R0/phi0
K <- Arec * SSBpR_out
h_out <- K/(4 + K)
R0_out <- (5*h_out-1)/(1-h_out)/Brec/SSBpR_out
} else {
Arec <- (5*h)^1.25/phi0
Brec <- 1.25 * log(5*h)/R0/phi0
h_out <- 0.2 * (Arec * SSBpR_out)^0.8
R0_out <- 1.25 * log(5 * h_out)/Brec/SSBpR_out
}
list(h = h_out, R0 = R0_out)
})
R0 <- new_SR$R0
h <- new_SR$h
}
SSB0 <- R0 * SSBpR_out
SSB <- apply(naa * exp(-spawn_time_frac * (Maa + faa)) * fecaa, c(1, 3), sum, na.rm = TRUE)
D <- SSB[,nyears]/SSB0
OM@Name <- Name
OM@M<-OM@L50<-OM@L50_95<-OM@L5<-OM@LFS<-OM@Vmaxlen <-c(1,1)
OM@a<-OM@b<-1
# Dimensions
OM@nsim<-nsim
OM@nyears<-nyears
OM@maxage<-maxage
OM@proyears<-proyears
OM@interval<-interval
OM <- simup(D, OM)
hs <- h
OM <- simup(hs, OM, do_slot = FALSE)
OM@h <- rep(hs[1], 2)
OM <- simup(R0, OM)
OM@Size_area_1 <- OM@Frac_area_1 <- OM@Prob_staying <- rep(0.5, 2)
# 'Contant terms'
LenCV=0.1
OM<-simup(LenCV,OM)
OM@SRrel <- SRrel # 1 = BevHolt, 2 = Ricker
OM@isRel<-FALSE # absolute selectivity relative to maturity - no used here
OM@CurrentYr<-CurrentYr
# Time varying terms
OM@Msd<-OM@Ksd<-OM@Linfsd<-OM@qcv<-OM@Esd<-OM@AC<-rep(0,2)
# Closed loop MP terms
OM@pstar<-0.5
OM@reps<-1
# dummy - filled to get past checks but not used due to cpars matrices in VPA
OM@EffYears<-c(1,5,OM@nyears)
OM@EffLower<-rep(1,3)
OM@EffUpper<-rep(1.1,3)
# Niche terms (assign defaults)
OM@Spat_targ<-c(1,1)
OM@qinc<-rep(0,2)
# Invent an OM with full observation error model for replacing
OM<-suppressMessages(Replace(OM, Obs, Sub = "Obs"))
OM<-suppressMessages(Replace(OM, Imp, Sub = "Imp"))
# Custom parameter arrays
Wt_age <- M_ageArray <- Len_age <- Mat_age <- V <- Fec_age <- array(NA, c(nsim, n_age, nyears + proyears))
# Historical filling
Find<-apply(faa,c(1,3),max)
Wt_age[,,1:nyears]<-waa
M_ageArray[,,1:nyears]<-Maa
Len_age[,,1:nyears]<-laa
Mat_age[,,1:nyears]<-Mataa
Fmax <-aperm(array(rep(Find,n_age),c(nsim,nyears,n_age)),c(1,3,2))
V[,,1:nyears]<-faa/Fmax
adjustV <- function(Vi) {
nan.ind <- apply(Vi, 2, max) %>% is.nan() %>% which()
if (length(nan.ind)>0) {
non.nan.ind <- max(nan.ind)+1
Vi[, nan.ind] <- Vi[, non.nan.ind]
}
Vi
}
V <- sapply(1:nsim, function(i) adjustV(V[i,,]), simplify = 'array') %>%
aperm(c(3,1,2))
# Future filling
parmu<-function(arr,nyears,proyears,nyr_par_mu){ # function for calculation mean values over last nyr_par_mu years
arr[,,nyears+(1:proyears)]<-array(rep(apply(arr[,,nyears-(0:(nyr_par_mu-1))],1:2,mean),proyears),c(dim(arr)[1:2],proyears))
arr
}
# Projected arrays based on mean of last nyr_par_mu years
Wt_age<-parmu(Wt_age,nyears,proyears,nyr_par_mu)
M_ageArray<-parmu(M_ageArray,nyears,proyears,nyr_par_mu)
Len_age<-parmu(Len_age,nyears,proyears,nyr_par_mu)
Mat_age<-parmu(Mat_age,nyears,proyears,nyr_par_mu)
V<-parmu(V,nyears,proyears,nyr_par_mu)
OM@cpars$Wt_age <- Wt_age
OM@cpars$M_ageArray <- M_ageArray
OM@cpars$Len_age <- Len_age
OM@cpars$Mat_age <- Mat_age
OM@cpars$V <- V
OM@cpars$Find <- Find
if(!is.null(dots$fecaa)) {
Fec_age[,,1:nyears] <- fecaa
OM@cpars$Fec_age <- parmu(Fec_age,nyears,proyears,nyr_par_mu)
}
# Deterministic Recruitment
if(SRrel == 1) {
recd <- vapply(1:nsim, function(x) (0.8*R0[x]*h[x]*SSB[x, ])/(0.2*SSB0[x]*(1-h[x]) + (h[x]-0.2)*SSB[x, ]),
numeric(nyears)) %>% t()
} else {
recd <- local({
a <- (5 * h)^1.25/SSBpR_out
b <- 1.25 * log(5 * h)/SSBpR_out/R0
vapply(1:nsim, function(x) a[x] * SSB[x, ] * exp(-b[x] * SSB[x, ]), numeric(nyears)) %>% t()
})
}
recdevs <- log(naa[, 1, ]/recd) # age zero
recdevs[is.na(recdevs)] <- 0
if(!is.null(dots$Perr)) {
if(length(dots$Perr) == 1) {
procsd <- rep(dots$Perr, nsim)
} else {
procsd <- dots$Perr
}
} else {
procsd <- apply(recdevs, 1, sd)
}
if(!is.null(dots$AC)) {
if(length(dots$AC) == 1) {
AC <- rep(dots$AC, nsim)
} else {
AC <- dots$AC
}
} else {
AC <- apply(recdevs, 1, function(x) stats::acf(x, plot = FALSE)$acf[2,1,1])
}
# Initial naa initialization options (VPAs apparently do this differently...)
surv <- lapply(1:nsim, calc_survival,
StockPars = list(M_ageArray = Maa, n_age = n_age),
plusgroup = plusgroup, inc_spawn_time = FALSE) %>%
abind(along = 3) %>%
aperm(c(3,1,2))
Perr <- array(NA_real_, c(nsim, maxage + nyears - LowerTri))
if(altinit < 2) { # normal assumption with or without plusgroup
Perr[, n_age:1] <- log(naa[, , 1]/(R0 * surv[, , 1]))
} else if(altinit==2) { # temporary fix for DLMtool initialization of plusgroup
Perr[, n_age:2] <- log(naa[, 1:(n_age-1), 1]/(R0 * SBsurv[, 1:(n_age-1), 1]))
survDLMtool<-aperm(exp(-apply(Maa[,,1],1,cumsum)),c(2,1))
fac<-surv[, n_age, 1]/survDLMtool[, n_age]
Perr[, 1] <- log(naa[, n_age, 1]/(R0 * surv[, n_age, 1] * fac))
}
Perr[, maxage + 2:(nyears - LowerTri)] <- recdevs[, 2:(nyears - LowerTri)]
Perr_pro <- sample_recruitment(Perr_hist = Perr, proyears = proyears + LowerTri, procsd = procsd, AC = AC)
OM@cpars$Perr_y <- exp(cbind(Perr, Perr_pro))
OM@Perr <- rep(mean(procsd),2)
OM@AC <- rep(mean(AC), 2)
OM@cpars$AC <- AC
OM@cpars$Perr <- procsd
if(fixq1) OM@cpars$qs <- rep(1, nsim) # Overrides q estimation to fix q at 1 for VPA for which F history is
OM@maxF <- ceiling(max(Find))
OM@Linf <- rep(1,2)
OM@K <- rep(0.3,2)
OM@t0 <- rep(0,2)
OM@DR <- rep(0,2)
OM@MPA<-FALSE
if (any(spawn_time_frac > 0)) OM@cpars$spawn_time_frac <- spawn_time_frac
if(!plusgroup) OM@cpars$plusgroup <- 0L
if(report){ # Produce a quick diagnostic plot of OM vs VPA numbers at age
if(!silent) message("\nRunning historical simulations to compare VPA output and OM conditioning...\n")
Hist <- runMSE(OM, Hist = TRUE, silent = TRUE)
nc<-ceiling((maxage+1)/3)
nr<-ceiling((maxage+1)/nc)
old_par <- par(no.readonly = TRUE)
on.exit(par(old_par))
par(mfrow=c(nr,nc),mai=c(0.4,0.4,0.3,0.05),omi=c(0.25,0.25,0.01,0.01))
yrs<-CurrentYr-((nyears-1):0)
cols<-rep('black',nyears)
pch<-rep(1,nyears)
cols[nyears-(0:(LowerTri-1))]<-"blue"
pch[nyears-(0:(LowerTri-1))]<-4
for(a in 0:maxage){
N_OM<-apply(Hist@AtAge$Number[1,a+1,,], 1, sum)
ylim=c(0,max(naa[1,a+1,],N_OM, na.rm = TRUE) * 1.05)
if(a==0) plot(yrs,naa[1,a+1,],xlab="",ylab="",col=cols,pch=pch,ylim=ylim,yaxs='i')
if(a>0) plot(yrs,naa[1,a+1,],xlab="",ylab="",ylim=ylim,yaxs='i')
lines(yrs,N_OM,col='green')
mtext(paste("Age ",a),3,line=0.5,cex=0.9)
if(a==0)legend('top',legend=c("Assessment","OM","Recr. ignored (LowerTri)"),text.col=c('black','green','blue'),bty='n')
res<-N_OM-naa[1,a+1,]
plotres<- abs(res) > mean(naa[1,a+1,]*0.025, na.rm = TRUE)
if(any(plotres, na.rm = TRUE)) {
for(y in 1:nyears) {
if(!is.na(plotres[y]) && plotres[y]) lines(rep(yrs[y], 2), c(naa[1,a+1,y],N_OM[y]),col='red')
}
}
}
mtext("Year",1,line=0.3,outer=T)
mtext("Numbers",2,line=0.4,outer=T)
}
return(OM)
}
#' @rdname Assess2OM
#' @export
VPA2OM <- Assess2OM
#' Reads bootstrap estimates from a stock assessment model into a multi-fleet operating model.
#'
#' @description A function that develops a multiple fleet operating model (\linkS4class{MOM}) and either models a unisex or 2-sex
#' stock from arrays of abundance, fishing mortality, and biological parameters. The user still
#' needs to parameterize most of the observation and implementation portions of the operating model.
#' @param Name Character string. The name of the multi-OM.
#' @param proyears Positive integer. The number of projection years for MSE.
#' @param interval Positive integer. The interval at which management procedures will update the management advice in \link[MSEtool]{multiMSE}, e.g., 1 = annual updates.
#' @param CurrentYr Positive integer. The current year (e.g., final year of fitting to data)
#' @param h The steepness of the stock-recruitment curve. Either a single numeric or a length nsim vector.
#' @param Obs Either a single observation model to be used for all sexes and populations (class \linkS4class{Obs}), or a list where
#' \code{Obs[[f]]} is the \linkS4class{Obs} object for fleet `f` (identical between sexes).
#' @param Imp Either a single implementation model to be used for all sexes and populations (class \linkS4class{Imp}), or a list where
#' \code{Imp[[f]]} is the \linkS4class{Obs} object for fleet `f` (identical between sexes).
#' @param naa Numbers-at-age by sex `[first age is age zero]`. Four-dimensional numeric array `[sim, ages, year, p]`. `[p]` indexes the population, where
#' `[p = 1]` for females and `[p = 2]` for males.
#' @param faa Fishing mortality rate-at-age by sex and fleet `[first age is age zero]`. Five-dimensional numeric array `[sim, ages, year, p, f]` where `[f]` indexes fishery fleet.
#' @param waa Weight-at-age `[first age is age zero]`. Four-dimensional numeric array `[sim, ages, year, p]`.
#' @param Mataa Maturity (spawning fraction)-at-age `[first age is age zero]`. Four-dimensional numeric array `[sim, ages, year, p]`.
#' @param Maa Natural mortality rate-at-age `[first age is age zero]`. Four-dimensional numeric array `[sim, ages, year, p]`.
#' @param laa Length-at-age `[first age is age zero]`. Four-dimensional numeric array `[sim, ages, year, p]`.
#' @param fecaa Fecundity at age `[first age is age zero]`. If missing, default fecundity is the product of maturity and weight at age.
#' @param nyr_par_mu Positive integer. The number of recent years that natural mortality, age vulnerability, weight, length and maturity parameters are averaged over for defining future projection conditions.
#' @param LowerTri Integer. The number of recent years for which model estimates of recruitment are ignored (not reliably estimated by the assessment)
#' @param recind Positive integer. The first age class that fish 'recruit to the fishery'. The default is 0 - ie the first position in the age dimension of naa is age zero
#' @param plusgroup Logical. Does the assessment assume that the oldest age class is a plusgroup?
#' @param altinit Integer. Various assumptions for how to set up the initial numbers. 0: standard, 1: no plus group, 2: temporary fix for MSEtool plus group initialization
#' @param fixq1 Logical. Should q be fixed (ie assume the F-at-age array faa is accurate?
#' @param report Logical, if TRUE, a diagnostic will be reported showing the matching of the OM reconstructed numbers at age vs the assessment.
#' @param silent Whether to silence messages to the console.
#' @param ... Additional arguments (for all, either a numeric or a length nsim vector):
#' \itemize{
#' \item \code{SRrel} Stock-recruit relationship. (\code{1} for Beverton-Holt (default), \code{2} for Ricker)
#' \item \code{R0} unfished recruitment
#' \item \code{phi0} unfished spawners per recruit associated with R0 and h. With time-varying parameters, openMSE uses the mean phi0
#' in the first \code{ageM} (age of 50 percent maturity) years for the stock-recruit relationship. \code{Assess2OM} will re-calculate R0 and h
#' in the operating model such that the stock-recruit \code{alpha} and \code{beta} parameters match values implied in the input.
#' \item \code{Perr} recruitment standard deviation (lognormal distribution) for sampling future recruitment
#' \item \code{AC} autocorrelation in future recruitment deviates.
#' }
#' @details Use a seed for the random number generator to sample future recruitment.
#' @return An object of class \linkS4class{MOM}.
#' @author Q. Huynh
#' @export
#' @seealso \link{SS2MOM} \link{multiMSE} \link{Assess2OM}
Assess2MOM <- function(Name = "MOM created by Assess2MOM",
proyears = 50, interval = 2, CurrentYr = as.numeric(format(Sys.Date(), "%Y")),
h = 0.999,
Obs = MSEtool::Imprecise_Unbiased, Imp = MSEtool::Perfect_Imp,
naa, faa, waa, Mataa, Maa, laa, fecaa,
nyr_par_mu = 3, LowerTri = 1,
recind = 0, plusgroup = TRUE, altinit = 0, fixq1 = TRUE,
report = FALSE, silent = FALSE, ...) {
np <- dim(naa)[4]
ny <- dim(naa)[3]
nf <- dim(faa)[5]
if (!silent) message("A ", np, "-sex MOM will be created with ", nf, " fishing fleets and ", ny, " historical years.")
if (inherits(Obs, "Obs")) Obs <- replicate(nf, Obs)
if (inherits(Imp, "Imp")) Imp <- replicate(nf, Imp)
if (missing(fecaa)) fecaa <- Mataa * waa
OMs <- lapply(1:np, function(p) {
lapply(1:nf, function(f) {
Assess2OM(Name="A fishery made by VPA2OM",
proyears = proyears, interval = interval, CurrentYr = CurrentYr,
h = h, Obs = Obs[[f]], Imp = Imp[[f]],
naa = naa[, , , p],
faa = faa[, , , p, f],
waa = waa[, , , p],
Mataa = Mataa[, , , p],
Maa = Maa[, , , p],
laa = laa[, , , p],
nyr_par_mu = nyr_par_mu,
LowerTri = LowerTri,
recind = recind,
plusgroup = plusgroup,
altinit = altinit,
fixq1 = fixq1,
report = FALSE, # Plots will be inaccurate anyway
silent = TRUE,
fecaa = fecaa[, , , p],
...)
})
})
MOM <- suppressMessages(new("MOM"))
slot_intersect <- intersect(slotNames("MOM"), slotNames("OM"))
for(i in slot_intersect) slot(MOM, i) <- slot(OMs[[1]][[1]], i)
MOM@Name <- Name
MOM@cpars <- list()
MOM@cpars <- lapply(1:np, function(p) {
lapply(1:nf, function(f) {
cp <- OMs[[p]][[f]]@cpars
cp$Perr_y <- OMs[[1]][[1]]@cpars$Perr_y
cp$D <- rep(1, length(cp$D))
return(cp)
})
})
MOM@Stocks <- lapply(1:np, function(p) {
ss <- SubOM(OMs[[p]][[1]], "Stock")
if (np == 1) {
ss@Name <- "Unisex"
} else if (p == 1) {
ss@Name <- "Female"
} else {
ss@Name <- "Male"
}
return(ss)
})
MOM@Fleets <- lapply(1:np, function(p) {
lapply(1:nf, function(f) {
ff <- SubOM(OMs[[p]][[f]], "Fleet")
ff@Name <- paste("Fleet", f)
return(ff)
})
})
MOM@Obs <- lapply(1:np, function(p) Obs)
MOM@Imps <- lapply(1:np, function(p) Imp)
MOM@CatchFrac <- lapply(1:np, function(p) {
Z <- apply(faa[, , ny, p, ], 1:2, sum) + Maa[, , ny, p] # nsim x n_age
CAA <- sapply(1:nf, function(f) naa[, , ny, p] * faa[, , ny, p, f] * (1 - exp(-Z))/Z,
simplify = "array") # nsim x n_age x nf
CB <- apply(CAA * replicate(nf, waa[, , ny, p]), c(1, 3), sum)
t(apply(CB, 1, function(x) x/sum(x)))
})
if (np == 2) {
if (!silent) message("Stock recruit parameters based on female biology and recruitment predicted by female SSB.")
MOM@SexPars <- list(SSBfrom = matrix(c(1, 0), 2, 2, byrow = TRUE))
MOM@Complexes <- list(c(1, 2))
}
return(MOM)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.