Nothing
#' read control file from SS version 3.24
#'
#' Read Stock Synthesis (version 3.24) control file into list object in R.
#' This function comes with its wrapper function SS_readctl
#' that calls SS_readctl_3.24 (this function) or SS_readctl_3.30
#' (to be available in future).
#'
#'
#' @template file
#' @param version Deprecated. SS version number. Currently only "3.24" or "3.30" are supported,
#' either as character or numeric values (noting that numeric 3.30 = 3.3).
#' @template readctl_vars
#' @param Nfleet number of fisheries in the model. This information is also not
#' explicitly available in control file
#' @param Nsurveys number of survey fleets in the model. This information is also not
#' explicitly available in control file
#' @param N_CPUE_obs numeric vector of length=Nfleet+Nsurveys containing number
#' of data points of each CPUE time series
#' @param use_datlist LOGICAL if TRUE, use datlist to derive parameters which
#' can not be determined from control file. Defaults to TRUE
#' @param datlist list or character. if list : produced from SS_writedat
#' or character : file name of dat file.
#' @param ptype include a column in the output indicating parameter type?
#' (Can be useful, but causes problems for SS_writectl.) Defaults to FALSE.
#' @author Yukio Takeuchi, Neil Klaer, Iago Mosqueira, and Kathryn Doering
#' @export
#' @seealso [SS_readctl()], [SS_readdat()]
#' [SS_readdat_3.24()],[SS_readdat_3.30()]
#' [SS_readstarter()], [SS_readforecast()],
#' [SS_writestarter()],
#' [SS_writeforecast()], [SS_writedat()]
SS_readctl_3.24 <- function(file,
verbose = FALSE,
echoall = lifecycle::deprecated(), # soft deprecate
version = lifecycle::deprecated(), # soft deprecate
use_datlist = TRUE,
datlist = "data.ss_new",
# Parameters that are not defined in control file:
nseas = NULL,
N_areas = NULL,
Nages = NULL,
Ngenders = lifecycle::deprecated(), # soft deprecate
Nsexes = NULL,
Npopbins = NA,
Nfleet = NULL,
Nsurveys = NULL,
Do_AgeKey = NULL,
N_tag_groups = NULL,
# This information is needed if Q_type of 3 or 4 is used
N_CPUE_obs = NULL,
##################################
ptype = FALSE) {
# deprecated variable warnings -----
# soft deprecated for now, but fully deprecate in the future.
if (lifecycle::is_present(echoall)) {
lifecycle::deprecate_warn(
when = "1.41.1",
what = "SS_readctl_3.24(echoall)"
)
}
if (lifecycle::is_present(version)) {
lifecycle::deprecate_warn(
when = "1.41.1",
what = "SS_readctl_3.24(version)"
)
}
version <- "3.24"
if (lifecycle::is_present(Ngenders)) {
lifecycle::deprecate_warn(
when = "1.41.1",
what = "SS_readctl_3.24(Ngenders)",
details = "Please use Nsexes instead. Ability to use Ngenders will be dropped in next release."
)
Nsexes <- Ngenders
}
# internally used fun definitions ----
# function to read Stock Synthesis data files
if (verbose) message("running SS_readctl_3.24\n")
dat <- readLines(file, warn = FALSE)
nver <- as.numeric(substring(version, 1, 4))
# parse all the numeric values into a long vector (allnums)
temp <- strsplit(dat[2], " ")[[1]][1]
if (!is.na(temp) && temp == "Start_time:") dat <- dat[-(1:2)]
allnums <- NULL
for (i in 1:length(dat)) {
# split along blank spaces
mysplit <- strsplit(dat[i], split = "[[:blank:]]+")[[1]]
mysplit <- mysplit[mysplit != ""]
# if final value is a number is followed immediately by a pound ("1#"),
# this needs to be split
nvals <- length(mysplit)
if (nvals > 0) mysplit[nvals] <- strsplit(mysplit[nvals], "#")[[1]][1]
# convert to numeric
nums <- suppressWarnings(as.numeric(mysplit))
if (sum(is.na(nums)) > 0) {
maxcol <- min((1:length(nums))[is.na(nums)]) - 1
} else {
maxcol <- length(nums)
}
if (maxcol > 0) {
nums <- nums[1:maxcol]
allnums <- c(allnums, nums)
}
}
# Function to add vector to ctllist
add_vec <- function(ctllist, length, name, comments = NULL) {
i <- ctllist$".i"
dat <- ctllist$".dat"
ctllist[["temp"]] <- dat[i + 1:length - 1]
ctllist$".i" <- i + length
if (is.null(comments)) {
names(ctllist[["temp"]]) <- paste0(
paste0("#_", name, "_", collapse = ""),
1:length
)
} else {
names(ctllist[["temp"]]) <- comments
}
if (!is.na(name)) names(ctllist)[names(ctllist) == "temp"] <- name
if (verbose) {
message(name, ", i = ", ctllist$".i", "\n")
print(ctllist[name])
}
return(ctllist)
}
# Function to add data as data.frame to ctllist
add_df <- function(ctllist, nrow, ncol, col.names, name, comments = NULL) {
k <- nrow * ncol
i <- ctllist$".i"
dat <- ctllist$".dat"
df0 <- as.data.frame(matrix(dat[i + 1:k - 1],
nrow = nrow, ncol = ncol,
byrow = TRUE
))
colnames(df0) <- col.names
if (is.null(comments)) {
rownames(df0) <- paste0(paste0("#_", name, collapse = ""), 1:nrow)
} else {
rownames(df0) <- comments
}
i <- i + k
ctllist[["temp"]] <- df0
ctllist$".i" <- i
if (!is.na(name)) names(ctllist)[names(ctllist) == "temp"] <- name
if (verbose) {
message(name, ", i = ", ctllist$".i", "\n")
print(ctllist[[which(names(ctllist) == name)]])
}
return(ctllist)
}
## function to add an element to ctllist
add_elem <- function(ctllist = NA, name) {
i <- ctllist$".i"
dat <- ctllist$".dat"
ctllist[["temp"]] <- dat[i]
ctllist$".i" <- i + 1
if (!is.na(name)) names(ctllist)[names(ctllist) == "temp"] <- name
if (verbose) {
message(
name, ", i = ", ctllist$".i", " ;",
ctllist[[which(names(ctllist) == name)]]
)
}
return(ctllist)
}
## function to add list to ctllist
add_list <- function(ctllist = NA, name, length, length_each) {
i <- ctllist$".i"
dat <- ctllist$".dat"
ctllist[["temp"]] <- list()
for (j in 1:length) {
ctllist[["temp"]][[j]] <- dat[i + 1:length_each[j] - 1]
i <- i + length_each[j]
}
ctllist$".i" <- i
if (!is.null(name)) names(ctllist)[names(ctllist) == "temp"] <- name
if (verbose) message(name, ", i = ", ctllist$".i")
return(ctllist)
}
# setup ----
# set initial position in the vector of numeric values
i <- 1
# create empty list to store quantities
ctllist <- list()
ctllist$".i" <- i
ctllist$".dat" <- allnums
ctllist[["warnings"]] <- ""
if (!use_datlist) {
ctllist[["nseas"]] <- nseas
ctllist[["N_areas"]] <- N_areas
ctllist[["Nages"]] <- Nages
ctllist[["Nsexes"]] <- Nsexes
ctllist[["Npopbins"]] <- Npopbins
ctllist[["Nfleet"]] <- Nfleet
ctllist[["Nsurveys"]] <- Nsurveys
ctllist[["Do_AgeKey"]] <- Do_AgeKey
ctllist[["N_tag_groups"]] <- N_tag_groups
ctllist[["N_CPUE_obs"]] <- N_CPUE_obs
fleetnames <- paste0("FL", 1:Nfleet)
if (Nsurveys > 0) fleetnames <- c(fleetnames, paste0("S", 1:Nsurveys))
ctllist[["fleetnames"]] <- fleetnames
} else {
if (is.character(datlist)) {
if (!file.exists(datlist)) {
stop("Cannot find data file specified in datlist: ", datlist)
}
datlist <- SS_readdat(file = datlist, version = "3.24", verbose = FALSE)
}
if (is.null(datlist)) {
stop("datlist from SS_readdat is needed if use_datlist is TRUE")
}
ctllist[["nseas"]] <- nseas <- datlist[["nseas"]]
ctllist[["N_areas"]] <- N_areas <- datlist[["N_areas"]]
ctllist[["Nages"]] <- Nages <- datlist[["Nages"]]
ctllist[["Nsexes"]] <- Nsexes <- datlist[["Ngenders"]] # note: change if renamed to datlist[["Nsexes"]]
ctllist[["Npopbins"]] <- Npopbins <- datlist[["N_lbinspop"]]
ctllist[["Nfleet"]] <- Nfleet <- datlist[["Nfleet"]]
ctllist[["Nsurveys"]] <- Nsurveys <- datlist[["Nsurveys"]]
if (datlist[["N_ageerror_definitions"]] > 0) {
ctllist[["Do_AgeKey"]] <- Do_AgeKey <-
ifelse(any(datlist[["ageerror"]][1:(nrow(datlist[["ageerror"]]) / 2) * 2, 1] < 0),
1,
0
)
}
ctllist[["N_tag_groups"]] <- N_tag_groups <- datlist[["N_tag_groups"]]
N_CPUE_obs <- lapply(
seq_len(Nfleet + Nsurveys),
function(i) {
sum(datlist[["CPUE"]][, "index"] == i)
}
)
N_CPUE_obs <- unlist(N_CPUE_obs)
ctllist[["N_CPUE_obs"]] <- N_CPUE_obs
ctllist[["fleetnames"]] <- fleetnames <- datlist[["fleetnames"]]
}
# specifications ----
ctllist[["sourcefile"]] <- file
ctllist[["type"]] <- "Stock_Synthesis_control_file"
ctllist[["ReadVersion"]] <- "3.24"
ctllist[["eof"]] <- FALSE
# internally used commmon values ----
lng_par_colnames <- c(
"LO", "HI", "INIT", "PRIOR", "PR_type", "SD", "PHASE",
"env_var", "use_dev", "dev_minyr", "dev_maxyr",
"dev_stddev", "Block", "Block_Fxn"
)
srt_par_colnames <- c("LO", "HI", "INIT", "PRIOR", "PR_type", "SD", "PHASE")
if (verbose) message("SS_readctl_3.24 - read version = ", ctllist[["ReadVersion"]])
# beginning of ctl ----
# model dimensions
ctllist <- add_elem(ctllist, "N_GP")
ctllist <- add_elem(ctllist, "N_platoon")
if (ctllist[["N_platoon"]] > 1) {
stop("sub morphs are not supported yet")
ctllist <- add_elem(ctllist, "submorphdist")
} else {
ctllist[["sd_ratio"]] <- 1.0
ctllist[["submorphdist"]] <- 1.0
}
if (ctllist[["submorphdist"]][1] < 0.0) {
if (ctllist[["N_platoon"]] == 1) {
ctllist[["submorphdist"]] <- 1.0
} else if (ctllist[["N_platoon"]] == 3) {
ctllist[["submorphdist"]] <- c(0.15, 0.70, 0.15)
} else if (ctllist[["N_platoon"]] == 5) {
ctllist[["submorphdist"]] <- c(0.031, 0.237, 0.464, 0.237, 0.031)
}
}
ctllist[["submorphdist"]] <- ctllist[["submorphdist"]] / sum(ctllist[["submorphdist"]])
if (ctllist[["N_GP"]] * ctllist[["nseas"]] * ctllist[["N_areas"]] > 1) {
ctllist <- add_elem(ctllist, "recr_dist_read")
recr_dist_read <- ctllist[["recr_dist_read"]]
# number of recruitment assignments (overrides GP*area*seas parameter values)
ctllist <- add_elem(ctllist, "recr_dist_inx")
# recruitment interaction requested
ctllist <- add_df(ctllist, "recr_dist_pattern",
nrow = recr_dist_read,
ncol = 3, col.names = c("GP", "seas", "area")
)
} else {
ctllist[["recr_dist_inx"]] <- FALSE
}
ctllist[["recr_dist_method"]] <- 1 # compatibility with v 3.30
if (ctllist[["N_areas"]] > 1) {
ctllist <- add_elem(ctllist, "N_moveDef") # _N_movement_definitions goes here if N_areas > 1
if (ctllist[["N_moveDef"]] > 0) {
ctllist <- add_elem(ctllist, "firstAgeMove") # _first age that moves (real age at begin of season, not integer) also cond on do_migration>0
ctllist <- add_df(ctllist, "moveDef",
nrow = ctllist[["N_moveDef"]], ncol = 6,
col.names = c("seas", "morph", "source", "dest", "age", "age2")
)
}
}
# block setup ----
ctllist <- add_elem(ctllist, "N_Block_Designs") # _Nblock_Patterns
if (ctllist[["N_Block_Designs"]] > 0) {
ctllist <- add_vec(ctllist,
name = "blocks_per_pattern",
length = ctllist[["N_Block_Designs"]]
)
# _blocks_per_pattern
# begin and end years of blocks
ctllist <- add_list(ctllist,
name = "Block_Design",
length = ctllist[["N_Block_Designs"]],
length_each = ctllist[["blocks_per_pattern"]] * 2
)
}
ctllist[["time_vary_adjust_method"]] <- 3 # compatibility with v 3.30
# compatibility with v 3.30
ctllist[["time_vary_auto_generation"]] <- c(1, 1, 1, 1, 1)
# MG setup ----
# M setup ----
ctllist <- add_elem(ctllist, "fracfemale") # _fracfemale
ctllist <- add_elem(ctllist, "natM_type") # _natM_type
if (ctllist[["natM_type"]] == 0) {
N_natMparms <- 1
} else if (ctllist[["natM_type"]] == 1) {
ctllist <- add_elem(ctllist, name = "N_natM") # _Number of M_segments
ctllist <- add_vec(ctllist,
name = "M_ageBreakPoints",
length = ctllist[["N_natM"]]
) # age(real) at M breakpoints
N_natMparms <- ctllist[["N_natM"]]
} else if (ctllist[["natM_type"]] == 2) {
N_natMparms <- 1
ctllist <- add_elem(ctllist, "Lorenzen_refage") # _natM_type
} else if (ctllist[["natM_type"]] %in% c(3, 4)) {
N_natMparms <- 0
ctllist <- add_df(ctllist,
name = "natM",
nrow = ctllist[["N_GP"]] * ctllist[["Nsexes"]],
ncol = Nages + 1,
col.names = paste0("Age_", 0:Nages)
)
} else {
stop("natM_type = ", ctllist[["natM_type"]], " is not yet implemented in this script")
}
if (verbose) message("N_natMparms = ", N_natMparms)
# growth setup ----
ctllist <- add_elem(ctllist, name = "GrowthModel")
# GrowthModel: 1=vonBert with L1&L2; 2=Richards with L1&L2; 3=age_specific_K; 4=not implemented
ctllist <- add_elem(ctllist, name = "Growth_Age_for_L1") # _Growth_Age_for_L1
ctllist <- add_elem(ctllist, name = "Growth_Age_for_L2") # _Growth_Age_for_L2 (999 to use as Linf)
ctllist[["Exp_Decay"]] <- 0.2 # Exponential decay for growth above maximum age - default 0.2
if (ctllist[["GrowthModel"]] == 3) {
ctllist <- add_elem(ctllist, name = "N_ageK")
}
if (ctllist[["GrowthModel"]] == 1) { # 1=vonBert with L1&L2
N_growparms <- 5
} else if (ctllist[["GrowthModel"]] == 2) { # 2=Richards with L1&L2
N_growparms <- 6
} else if (ctllist[["GrowthModel"]] == 3) { # 3=age_specific_K
N_growparms <- 5 + ctllist[["N_ageK"]]
ctllist <- add_vec(ctllist, name = "Age_K_points", length = ctllist[["N_ageK"]])
# points at which age-specific multipliers to K will be applied
} else if (ctllist[["GrowthModel"]] == 4) {
stop("GrowthModel == 4 is not implemented")
N_growparms <- 2 # for the two CV parameters
k1 <- ctllist[["N_GP"]] * ctllist[["Nsexes"]] # for reading age_natmort
ctllist <- add_df(ctllist,
name = "Len_At_Age_rd", nrow = k1,
ncol = Nages + 1, col.names = paste0("Age_", 0:Nages)
)
} else {
stop("GrowthModel ", ctllist[["GrowthModel"]], " is not supported yet.")
}
MGparm_per_def <- N_natMparms + N_growparms
ctllist[["N_natMparms"]] <- N_natMparms
ctllist <- add_elem(ctllist, name = "SD_add_to_LAA") # _SD_add_to_LAA (set to 0.1 for SS2 V1.x compatibility)
ctllist <- add_elem(ctllist, name = "CV_Growth_Pattern")
# _CV_Growth_Pattern: 0 CV=f(LAA); 1 CV=F(A); 2 SD=F(LAA); 3 SD=F(A); 4 logSD=F(A)
# maturity and MG options setup ----
ctllist <- add_elem(ctllist, name = "maturity_option")
# maturity_option: 1=length logistic; 2=age logistic;
# 3=read age-maturity by GP; 4=read age-fecundity by GP; 5=read fec and wt
# from wtatage.ss; 6=read length-maturity by GP
if (ctllist[["maturity_option"]] == 6) {
ctllist[["EmpiricalWAA"]] <- 1
} else {
ctllist[["EmpiricalWAA"]] <- 0
}
if (ctllist[["maturity_option"]] %in% c(3, 4)) {
ctllist <- add_df(ctllist,
name = "Age_Maturity", nrow = ctllist[["N_GP"]],
ncol = Nages + 1, col.names = paste0("Age_", 0:Nages)
)
} else if (ctllist[["maturity_option"]] == 6) {
ctllist <- add_df(ctllist,
name = "Length_Maturity", nrow = ctllist[["N_GP"]],
ncol = Npopbins, col.names = paste0("Len_", 1:Npopbins)
)
}
ctllist <- add_elem(ctllist, "First_Mature_Age")
ctllist <- add_elem(ctllist, "fecundity_option")
ctllist <- add_elem(ctllist, "hermaphroditism_option")
if (ctllist[["hermaphroditism_option"]] > 0) {
ctllist <- add_elem(ctllist, "Herm_season")
ctllist <- add_elem(ctllist, "Herm_MalesInSSB")
}
ctllist <- add_elem(ctllist, "parameter_offset_approach")
ctllist <- add_elem(ctllist, "env_block_dev_adjust_method")
# MG parlines -----
N_MGparm <- MGparm_per_def * ctllist[["N_GP"]] * ctllist[["Nsexes"]]
MGparmLabel <- list()
cnt <- 1
# store parameter types M=1, Growth=2, WtLn = 3, Maturity = 4, Fecundity = 5,
# Hermaph = 6, RecDevs GP = 7 Areas = 8 Seas= 9, RecDev Interactions = 10,
# GrowthDevs = 11, Movement = 12
PType <- array()
GenderLabel <- c("Fem", "Mal")
for (i in 1:ctllist[["Nsexes"]]) {
for (j in 1:ctllist[["N_GP"]]) {
if (N_natMparms > 0) {
MGparmLabel[1:N_natMparms + cnt - 1] <-
paste0("NatM_p_", 1:N_natMparms, "_", GenderLabel[i], "_GP_", j)
PType[cnt:(N_natMparms + cnt - 1)] <- 1
cnt <- cnt + N_natMparms
}
if (ctllist[["GrowthModel"]] == 1) { # VB
tmp <- c(
"L_at_Amin_", "L_at_Amax_", "VonBert_K_", "CV_young_",
"CV_old_"
)
MGparmLabel[1:5 + cnt - 1] <- paste0(tmp, "_", GenderLabel[i], "_GP_", j)
PType[cnt:(5 + cnt - 1)] <- 2
cnt <- cnt + 5
} else if (ctllist[["GrowthModel"]] == 2) { # Richards
tmp <- c(
"L_at_Amin_", "L_at_Amax_", "VonBert_K_", "Richards_", "CV_young_",
"CV_old_"
)
MGparmLabel[1:6 + cnt - 1] <- paste0(tmp, "_", GenderLabel[i], "_GP_", j)
PType[cnt:(6 + cnt - 1)] <- 2
cnt <- cnt + 6
} else if (ctllist[["GrowthModel"]] == 3) {
tmp <- c(
"L_at_Amin_", "L_at_Amax_", "VonBert_K_",
paste0("Age_K_", 1:ctllist[["N_ageK"]]), "CV_young_", "CV_old_"
)
MGparmLabel[1:(5 + ctllist[["N_ageK"]]) + cnt - 1] <-
paste0(tmp, "_", GenderLabel[i], "_GP_", j)
PType[cnt:((5 + ctllist[["N_ageK"]]) + cnt - 1)] <- 2
cnt <- cnt + 5 + ctllist[["N_ageK"]]
}
}
}
N_MGparm <- N_MGparm + 2 * ctllist[["Nsexes"]] + 2 + 2 # add for wt-len(by gender), mat-len parms; eggs
MGparmLabel[cnt] <- paste0("Wtlen_1_", GenderLabel[1])
PType[cnt] <- 3
cnt <- cnt + 1
MGparmLabel[cnt] <- paste0("Wtlen_2_", GenderLabel[1])
PType[cnt] <- 3
cnt <- cnt + 1
MGparmLabel[cnt] <- paste0("Mat50%_", GenderLabel[1])
PType[cnt] <- 4
cnt <- cnt + 1
MGparmLabel[cnt] <- paste0("Mat_slope_", GenderLabel[1])
PType[cnt] <- 4
cnt <- cnt + 1
if (ctllist[["fecundity_option"]] == 1) {
MGparmLabel[cnt] <- paste0("Eggs/kg_inter_", GenderLabel[1])
PType[cnt] <- 5
cnt <- cnt + 1
MGparmLabel[cnt] <- paste0("Eggs/kg_slope_wt_", GenderLabel[1])
PType[cnt] <- 5
cnt <- cnt + 1
} else if (ctllist[["fecundity_option"]] == 2) {
MGparmLabel[cnt] <- paste0("Eggs_scalar_", GenderLabel[1])
PType[cnt] <- 5
cnt <- cnt + 1
MGparmLabel[cnt] <- paste0("Eggs_exp_len_", GenderLabel[1])
PType[cnt] <- 5
cnt <- cnt + 1
} else if (ctllist[["fecundity_option"]] == 3) {
MGparmLabel[cnt] <- paste0("Eggs_scalar_", GenderLabel[1])
PType[cnt] <- 5
cnt <- cnt + 1
MGparmLabel[cnt] <- paste0("Eggs_exp_wt_", GenderLabel[1])
PType[cnt] <- 5
cnt <- cnt + 1
} else if (ctllist[["fecundity_option"]] == 4) {
MGparmLabel[cnt] <- paste0("Eggs_intercept_", GenderLabel[1])
PType[cnt] <- 5
cnt <- cnt + 1
MGparmLabel[cnt] <- paste0("Eggs_slope_len_", GenderLabel[1])
PType[cnt] <- 5
cnt <- cnt + 1
} else if (ctllist[["fecundity_option"]] == 5) {
MGparmLabel[cnt] <- paste0("Eggs_intercept_", GenderLabel[1])
PType[cnt] <- 5
cnt <- cnt + 1
MGparmLabel[cnt] <- paste0("Eggs_slope_wt_", GenderLabel[1])
PType[cnt] <- 5
cnt <- cnt + 1
} else if (ctllist[["fecundity_option"]] == 6) { # check what to do with option 6
MGparmLabel[cnt] <- paste0("Eggs_intercept_", GenderLabel[1])
PType[cnt] <- 5
cnt <- cnt + 1
MGparmLabel[cnt] <- paste0("Eggs_slope_wt_", GenderLabel[1])
PType[cnt] <- 5
cnt <- cnt + 1
} else {
stop("Fecundity option ", ctllist[["fecundity_option"]], " is not supported")
}
if (ctllist[["Nsexes"]] == 2) {
MGparmLabel[cnt] <- paste0("Wtlen_1_", GenderLabel[2])
PType[cnt] <- 3
cnt <- cnt + 1
MGparmLabel[cnt] <- paste0("Wtlen_2_", GenderLabel[2])
PType[cnt] <- 3
cnt <- cnt + 1
}
if (ctllist[["hermaphroditism_option"]]) {
MGparmLabel[cnt] <- paste0("Herm_Infl_age", GenderLabel[1])
PType[cnt] <- 6
cnt <- cnt + 1
MGparmLabel[cnt] <- paste0("Herm_stdev", GenderLabel[1])
PType[cnt] <- 6
cnt <- cnt + 1
MGparmLabel[cnt] <- paste0("Herm_asymptote", GenderLabel[1])
PType[cnt] <- 6
cnt <- cnt + 1
N_MGparm <- N_MGparm + 3
}
N_MGparm <- N_MGparm + ctllist[["N_GP"]] + N_areas + nseas # add for the assignment to areas
MGparmLabel[cnt + 1:ctllist[["N_GP"]] - 1] <- paste0("RecrDist_GP_", 1:ctllist[["N_GP"]])
PType[cnt:(cnt + ctllist[["N_GP"]] - 1)] <- 7
cnt <- cnt + ctllist[["N_GP"]]
MGparmLabel[cnt + 1:N_areas - 1] <- paste0("RecrDist_Area_", 1:N_areas)
PType[cnt:(cnt + N_areas)] <- 8
cnt <- cnt + N_areas
MGparmLabel[cnt + 1:nseas - 1] <- paste0("RecrDist_Seas_", 1:nseas)
PType[cnt:(cnt + nseas - 1)] <- 9
cnt <- cnt + nseas
## number of recruiment distribution parameters
N_RecrDist_parms <- ctllist[["N_GP"]] + ctllist[["N_areas"]] + ctllist[["nseas"]]
if (ctllist[["recr_dist_inx"]]) {
## Interactions
N_MGparm <- N_MGparm + ctllist[["N_GP"]] * N_areas * nseas
N_RecrDist_parms <- N_RecrDist_parms + ctllist[["N_GP"]] * N_areas * nseas
for (i in 1:ctllist[["N_GP"]]) {
for (j in 1:nseas) {
for (k in 1:N_areas) {
MGparmLabel[cnt] <-
paste0("RecrDist_interaction_GP_", i, "_seas_", j, "_area_", k)
PType[cnt] <- 10
cnt <- cnt + 1
}
}
}
# add for the morph assignments within each area
}
# MGparms
N_MGparm <- N_MGparm + 1 # add 1 parameter for cohort-specific growth parameter
MGparmLabel[cnt] <- "CohortGrowDev"
PType[cnt] <- 11
cnt <- cnt + 1
if ((N_areas > 1) && (ctllist[["N_moveDef"]] > 0)) {
N_MGparm <- N_MGparm + ctllist[["N_moveDef"]] * 2
for (i in 1:ctllist[["N_moveDef"]]) {
seas <- ctllist[["moveDef"]][i, 1]
GP <- ctllist[["moveDef"]][i, 2]
from <- ctllist[["moveDef"]][i, 3]
to <- ctllist[["moveDef"]][i, 4]
MGparmLabel[cnt + 0:1] <- paste0(
"MoveParm_", c("A", "B"), "_seas_", seas,
"_GP_", GP, "_from_", from, "_to_", to
)
PType[cnt:(cnt + 1)] <- 12
cnt <- cnt + 2
}
}
if (Do_AgeKey) {
MGparmLabel[cnt + 0:6] <- paste0("AgeKeyParm", 1:7)
PType[cnt:(cnt + 6)] <- 13
cnt <- cnt + 7
N_MGparm <- N_MGparm + 7
}
ctllist <- add_df(ctllist,
name = "MG_parms", nrow = N_MGparm, ncol = 14,
col.names = lng_par_colnames, comments = MGparmLabel
)
ctllist[["MG_parms"]] <- cbind(ctllist[["MG_parms"]], PType)
# MG timevarying parlines ------
# environmental linkage lines
if (any(ctllist[["MG_parms"]][["env_var"]] != 0)) {
ctllist <- add_elem(ctllist, "read_MG_custom_env_var")
if (ctllist[["read_MG_custom_env_var"]] == 1) {
tmp_nrow <- length(ctllist[["MG_parms"]][["env_var"]][ctllist[["MG_parms"]][["env_var"]] != 0])
tmp_lab <- paste0("MG_env_var_", seq_len(tmp_nrow))
ctllist <- add_df(ctllist,
name = "MG_custom_env_var", nrow = tmp_nrow,
ncol = length(srt_par_colnames),
col.names = srt_par_colnames,
comments = tmp_lab
)
}
}
# time block parameters
nbp <- 0
cnt <- 1
PType <- array()
MGblockLabel <- list()
for (i in 1:N_MGparm) {
if (ctllist[["MG_parms"]][i, ][["Block"]] > 0) {
nv <- as.numeric(ctllist[["blocks_per_pattern"]][ctllist[["MG_parms"]][i, ][["Block"]]])
MGblockLabel[(nbp + 1):(nbp + nv)] <-
paste0(rownames(ctllist[["MG_parms"]][i, ]), "#", (nbp + 1):(nbp + nv))
nbp <- nbp + nv
PType[cnt:(nv + cnt - 1)] <- 15
cnt <- cnt + nv
}
}
if (nbp > 0) {
ctllist <- add_elem(ctllist, "custom_MG_block")
ctllist <- add_df(ctllist,
name = "MG_parms_blocks", nrow = nbp, ncol = 7,
col.names = srt_par_colnames, comments = MGblockLabel
)
ctllist[["MG_parms_blocks"]] <- cbind(ctllist[["MG_parms_blocks"]], PType)
}
# Read seasonal effects ----
ctllist <- add_vec(ctllist, name = "MGparm_seas_effects", length = 10)
PType <- array()
N_seas_effects <-
length(ctllist[["MGparm_seas_effects"]][ctllist[["MGparm_seas_effects"]] > 0])
if (N_seas_effects > 0) {
ctllist <- add_df(ctllist, "MG_parms_seas",
nrow = N_seas_effects * ctllist[["nseas"]], ncol = 7,
col.names = srt_par_colnames
)
PType[1:(N_seas_effects * ctllist[["nseas"]])] <- 16
ctllist[["MG_parms_seas"]] <- cbind(ctllist[["MG_parms_seas"]], PType)
}
DoParmDev <- sum(ctllist[["MG_parms"]][, 9])
if (DoParmDev > 0) {
ctllist <- add_elem(ctllist, "MGparm_Dev_Phase")
}
# SR ----
ctllist <- add_elem(ctllist, "SR_function")
N_SRparm <- c(0, 2, 2, 2, 3, 2, 3, 3, 0, 0)
N_SRparm2 <- N_SRparm[as.numeric(ctllist[["SR_function"]])] + 4
if (is.na(ctllist[["SR_function"]])) {
stop("SR_function is NA")
}
# SR parms ----
SRparmsLabels <- if (ctllist[["SR_function"]] == 3) {
# B-H SRR
c(
"SR_LN(R0)", "SR_BH_steep", "SR_sigmaR", "SR_envlink", "SR_R1_offset",
"SR_autocorr"
)
} else if (ctllist[["SR_function"]] == 2) {
# Ricker SRR
c(
"SR_LN(R0)", "SR_Ricker", "SR_sigmaR", "SR_envlink", "SR_R1_offset",
"SR_autocorr"
) ## Need to rivise with example inputs
} else if (ctllist[["SR_function"]] == 4) {
# SCAA
c(
"SR_LN(R0)", "SR_SCAA_null", "SR_sigmaR", "SR_envlink", "SR_R1_offset",
"SR_autocorr"
)
} else if (ctllist[["SR_function"]] == 5) {
# Hockey stick
c(
"SR_LN(R0)", "SR_hockey_infl", "SR_hockey_min_R", "SR_sigmaR",
"SR_envlink", "SR_R1_offset", "SR_autocorr"
)
} else if (ctllist[["SR_function"]] == 6) {
# B-H-flat SRR
c(
"SR_LN(R0)", "SR_BH_flat_steep", "SR_sigmaR", "SR_envlink",
"SR_R1_offset", "SR_autocorr"
)
} else if (ctllist[["SR_function"]] == 7) {
# survival_3Parm
c(
"SR_LN(R0)", "SR_surv_Sfrac", "SR_surv_Beta", "SR_sigmaR", "SR_envlink",
"SR_R1_offset", "SR_autocorr"
)
} else if (ctllist[["SR_function"]] == 8) {
# Shepard_3Parm
c(
"SR_LN(R0)", "SR_steepness", "SR_Shepard_c", "SR_sigmaR", "SR_envlink",
"SR_R1_offset", "SR_autocorr"
)
} else {
stop("SR_function ", ctllist[["SR_function"]], " is not supported yet.")
}
PType <- array()
ctllist <- add_df(ctllist,
name = "SR_parms", nrow = N_SRparm2, ncol = 7,
col.names = srt_par_colnames, comments = SRparmsLabels
)
PType[1:N_SRparm2] <- 17
ctllist[["SR_parms"]] <- cbind(ctllist[["SR_parms"]], PType)
# SR timevarying parlines ----
ctllist <- add_elem(ctllist, "SR_env_link") # _SR_env_link
ctllist <- add_elem(ctllist, "SR_env_target") # _SR_env_target_0=none;1=devs;_2=R0;_3=steepness
# recdevs ----
ctllist <- add_elem(ctllist, "do_recdev") # do_recdev: 0=none; 1=devvector; 2=simple deviations
ctllist <- add_elem(ctllist, "MainRdevYrFirst") # first year of main recr_devs; early devs can preceed this era
ctllist <- add_elem(ctllist, "MainRdevYrLast") # last year of main recr_devs; forecast devs start in following year
ctllist <- add_elem(ctllist, "recdev_phase") # _recdev phase
ctllist <- add_elem(ctllist, "recdev_adv") # (0/1) to read 13 advanced options
if (ctllist[["recdev_adv"]]) {
ctllist <- add_elem(ctllist, "recdev_early_start") # _recdev_early_start (0=none; neg value makes relative to recdev_start)
ctllist <- add_elem(ctllist, "recdev_early_phase") # _recdev_early_phase
ctllist <- add_elem(ctllist, "Fcast_recr_phase") # _forecast_recruitment phase (incl. late recr) (0 value resets to maxphase+1)
ctllist <- add_elem(ctllist, "lambda4Fcast_recr_like") # _lambda for Fcast_recr_like occurring before endyr+1
ctllist <- add_elem(ctllist, "last_early_yr_nobias_adj") # _last_early_yr_nobias_adj_in_MPD
ctllist <- add_elem(ctllist, "first_yr_fullbias_adj") # _first_yr_fullbias_adj_in_MPD
ctllist <- add_elem(ctllist, "last_yr_fullbias_adj") # _last_yr_fullbias_adj_in_MPD
ctllist <- add_elem(ctllist, "first_recent_yr_nobias_adj") # _first_recent_yr_nobias_adj_in_MPD
ctllist <- add_elem(ctllist, "max_bias_adj") # _max_bias_adj_in_MPD (-1 to override ramp and set biasadj=1.0 for all estimated recdevs)
ctllist <- add_elem(ctllist, "period_of_cycles_in_recr") # _period of cycles in recruitment (N parms read below)
ctllist <- add_elem(ctllist, "min_rec_dev") # min rec_dev
ctllist <- add_elem(ctllist, "max_rec_dev") # max rec_dev
ctllist <- add_elem(ctllist, "N_Read_recdevs") # _read_recdevs
if (ctllist[["period_of_cycles_in_recr"]] > 0) {
stop("Reading full parameters for recr cycles is not yet coded")
}
if (ctllist[["N_Read_recdevs"]] > 0) {
ctllist <- add_df(ctllist, "recdev_input",
ncol = 2,
nrow = ctllist[["N_Read_recdevs"]],
col.names = c("Year", "recdev")
)
}
}
# F setup ----
ctllist <- add_elem(ctllist, "F_ballpark") # F ballpark for annual F (=Z-M) for specified year
ctllist <- add_elem(ctllist, "F_ballpark_year") # F ballpark year (neg value to disable)
ctllist <- add_elem(ctllist, "F_Method") # F_Method: 1=Pope; 2=instan. F; 3=hybrid (hybrid is recommended)
ctllist <- add_elem(ctllist, "maxF") # max F or harvest rate, depends on F_Method
if (ctllist[["F_Method"]] == 1) {
# stop("stop currently F_method:1 is not implemented")
} else if (ctllist[["F_Method"]] == 2) {
# stop("stop currently F_method:2 is not implemented")
ctllist <- add_vec(ctllist, "F_setup", length = 3) # overall start F value; overall phase; N detailed inputs to read
if (ctllist[["F_setup"]][3] > 0) {
ctllist <- add_df(ctllist,
name = "F_setup2", nrow = ctllist[["F_setup"]][3],
ncol = 6,
col.names = c(
"fleet", "yr", "seas", "Fvalue", "se",
"phase"
)
)
}
} else if (ctllist[["F_Method"]] == 3) {
ctllist <- add_elem(ctllist, "F_iter") # N iterations for tuning F in hybrid method (recommend 3 to 7)
}
#
# _initial_F_parms
comments_initF <- paste0("InitF_", 1:Nfleet, "_", fleetnames[1:Nfleet])
PType <- array()
ctllist <- add_df(ctllist,
name = "init_F", nrow = Nfleet, ncol = 7,
col.names = srt_par_colnames, comments = comments_initF
)
PType[1:Nfleet] <- 18
ctllist[["init_F"]] <- cbind(ctllist[["init_F"]], PType)
# Q_setup ----
comments_fl <- paste0(1:(Nfleet + Nsurveys), " ", fleetnames)
ctllist <- add_df(ctllist,
name = "Q_setup", nrow = Nfleet + Nsurveys,
ncol = 4,
col.names = c("Den_dep", "env_var", "extra_se", "Q_type"),
comments = comments_fl
)
Do_Q_detail <- FALSE
if (sum(ctllist[["Q_setup"]][, 4] %in% c(3, 4)) > 0) {
ctllist <- add_elem(ctllist, name = "Do_Q_detail") ##
Do_Q_detail <- ctllist[["Do_Q_detail"]]
}
nqp <- 0
comments_Q_type <- list()
# Density dependant Q(Q-power)
for (i in 1:(Nfleet + Nsurveys)) {
if (ctllist[["Q_setup"]][i, ][["Den_dep"]] > 0) {
comments_Q_type[(nqp + 1)] <-
paste0("Q_power_", i, "_", fleetnames[i], collapse = "")
nqp <- nqp + 1
}
}
# Q-env
for (i in 1:(Nfleet + Nsurveys)) {
if (ctllist[["Q_setup"]][i, ][["env_var"]] > 0) {
comments_Q_type[(nqp + 1)] <- paste0("Q_env_", i, "_", fleetnames[i],
collapse = ""
)
nqp <- nqp + 1
}
}
# Q_extraSD
for (i in 1:(Nfleet + Nsurveys)) {
if (ctllist[["Q_setup"]][i, ][["extra_se"]] > 0) {
comments_Q_type[(nqp + 1)] <- paste0("Q_extraSD_", i, "_", fleetnames[i],
collapse = ""
)
nqp <- nqp + 1
}
}
# q parlines ----
# Q_type options: <0=mirror, 0=float_nobiasadj, 1=float_biasadj, 2=parm_nobiasadj, 3=parm_w_random_dev, 4=parm_w_randwalk, 5=mean_unbiased_float_assign_to_parm
## currently only float_nobiasadj (Q_type==0) is suppoerted
# _for_env-var:_enter_index_of_the_env-var_to_be_linked
# _Den-dep env-var extra_se Q_type
## Then check if random Q parameters are used.
## If yes, read 1 number for flag to see if to read single parameter for each random Q or
## one parameter for each data point
# Q-type
N_Q_parms <- nqp
i <- nqp + 1
for (fl in 1:(Nfleet + Nsurveys)) {
flname <- fleetnames[fl]
.Q_type <- ctllist[["Q_setup"]][fl, 4]
if (.Q_type == 2) { # One Q parameter
N_Q_parms <- N_Q_parms + 1
comments_Q_type[[i]] <-
paste0("LnQ_base_", fl, "_", flname, collapse = "")
i <- i + 1
} else if (.Q_type == 3) { # Random Q deviations
N_Q_parms <- N_Q_parms + 1
comments_Q_type[[i]] <- paste0("LnQ_base_", fl, "_", flname,
collapse = ""
)
i <- i + 1
if (Do_Q_detail) {
for (j in 1:(N_CPUE_obs[fl])) {
comments_Q_type[[i + j - 1]] <- paste0("Q_dev_", j, "_", fl, "_", flname)
}
N_Q_parms <- N_Q_parms + N_CPUE_obs[fl]
i <- i + j
}
} else if (.Q_type == 4) { # Random walk W
N_Q_parms <- N_Q_parms + 1
comments_Q_type[[i]] <- paste0("LnQ_base_", fl, "_", flname,
collapse = ""
)
i <- i + 1
if (Do_Q_detail) {
for (j in 1:(N_CPUE_obs[fl] - 1)) {
comments_Q_type[[i + j - 1]] <- paste0("Q_walk_", j, "_", fl, "_", flname)
}
N_Q_parms <- N_Q_parms + N_CPUE_obs[fl] - 1
i <- i + j
}
} else if (.Q_type == 5) {
N_Q_parms <- N_Q_parms + 1
comments_Q_type[[i]] <- paste0("LnQ_base_", fl, "_", flname,
collapse = ""
)
i <- i + 1
}
}
if (N_Q_parms > 0) {
ctllist <- add_df(ctllist,
name = "Q_parms", nrow = N_Q_parms, ncol = 7,
col.names = srt_par_colnames,
comments = unlist(comments_Q_type)
)
}
# selecitivty -----
# size_selex_types
comments_selex_types <- fleetnames
ctllist <- add_df(ctllist,
name = "size_selex_types",
nrow = Nfleet + Nsurveys, ncol = 4,
col.names = c("Pattern", "Discard", "Male", "Special"),
comments = comments_selex_types
)
size_selex_pattern_vec <- as.vector(ctllist[["size_selex_types"]][, 1])
ctllist <- add_df(ctllist,
name = "age_selex_types", nrow = Nfleet + Nsurveys,
ncol = 4,
col.names = c("Pattern", "Discard", "Male", "Special"),
comments = comments_selex_types
)
age_selex_pattern_vec <- as.vector(ctllist[["age_selex_types"]][, 1])
# 0 1 2 3 4 5 6 7 8 9 10
selex_patterns <- c(
0, 2, 8, 6, 0, 2, 2, 8, 8, 6, 0,
# 11 12 13 14 15 16 17 18 19 20
2, 2, 8, Nages + 1, 0, 2, Nages + 1, 8, 6, 6,
# 21 22 23 24 25 26 27 28 29 30
NA, 4, 6, 6, 3, 3, 3, NA, NA, 0,
# 31 32 33 34
0, 0, 0, 0
)
#########################################################
size_selex_Nparms <- vector(mode = "numeric", length = Nfleet + Nsurveys)
size_selex_label <- list()
for (j in 1:(Nfleet + Nsurveys)) {
# size_selex_pattern_vec
size_selex_Nparms[j] <- selex_patterns[size_selex_pattern_vec[j] + 1]
## some patterns need special treatment
special_selex <- FALSE # change flag to TRUE if need special treatment
if (size_selex_pattern_vec[j] == 27) { # spline
special_selex <- TRUE
size_selex_Nparms[j] <- size_selex_Nparms[j] +
ctllist[["size_selex_types"]][j, 4] * 2
size_selex_label[[j]] <- c(
paste0("SizeSpline_Code_", fleetnames[j], "_", j),
paste0("SizeSpline_GradLo_", fleetnames[j], "_", j),
paste0("SizeSpline_GradHi_", fleetnames[j], "_", j),
paste0("SizeSpline_Knot_", 1:ctllist[["size_selex_types"]][j, 4], "_", fleetnames[j], "_", j),
paste0("SizeSpline_Val_", 1:ctllist[["size_selex_types"]][j, 4], "_", fleetnames[j], "_", j)
)
}
if (!special_selex) {
size_selex_label[[j]] <- if (size_selex_Nparms[j] > 0) {
paste0("SizeSel_", j, "P_", 1:size_selex_Nparms[j], "_", fleetnames[j])
} else {
NULL
}
}
# do extra retention parameters
if ((ctllist[["size_selex_types"]][j, 2] > 0) && (ctllist[["size_selex_types"]][j, 2] <= 2)) # has discard type 1 or 2
{
if (size_selex_Nparms[j] > 0) {
size_selex_label[[j]] <- c(size_selex_label[[j]], paste0("SizeSel_", j, "PRet_", 1:4, "_", fleetnames[j]))
} else {
size_selex_label[[j]] <- paste0("SizeSel_", j, "PRet_", 1:4, "_", fleetnames[j])
}
size_selex_Nparms[j] <- size_selex_Nparms[j] + 4
}
if (ctllist[["size_selex_types"]][j, 2] == 2) # has discard type 2
{
if (size_selex_Nparms[j] > 0) {
size_selex_label[[j]] <- c(size_selex_label[[j]], paste0("SizeSel_", j, "PDis_", 1:4, "_", fleetnames[j]))
} else {
size_selex_label[[j]] <- paste0("SizeSel_", j, "PDis_", 1:4, "_", fleetnames[j])
}
size_selex_Nparms[j] <- size_selex_Nparms[j] + 4
}
# no extra retention parameters for discard type 3
# do extra offset parameters
if (ctllist[["size_selex_types"]][j, 3] == 1) # has value 1
{
if (size_selex_Nparms[j] > 0) {
size_selex_label[[j]] <- c(size_selex_label[[j]], paste0("SizeSel_", j, "PMale_", 1:4, "_", fleetnames[j]))
} else {
size_selex_label[[j]] <- paste0("SizeSel_", j, "PMale_", 1:4, "_", fleetnames[j])
}
size_selex_Nparms[j] <- size_selex_Nparms[j] + 4
}
if (ctllist[["size_selex_types"]][j, 3] == 2) # has value 2
{
if (size_selex_Nparms[j] > 0) {
size_selex_label[[j]] <- c(size_selex_label[[j]], paste0("SizeSel_", j, "PFemOff_", 1:4, "_", fleetnames[j]))
} else {
size_selex_label[[j]] <- paste0("SizeSel_", j, "PFemOff_", 1:4, "_", fleetnames[j])
}
size_selex_Nparms[j] <- size_selex_Nparms[j] + 4
}
if (ctllist[["size_selex_types"]][j, 3] == 3) # has value 3 - differs by version
{
if (nver < 3.24) {
nparms <- 3
} else {
nparms <- 5
}
if (size_selex_Nparms[j] > 0) {
size_selex_label[[j]] <- c(size_selex_label[[j]], paste0("SizeSel_", j, "PMalOff_", 1:nparms, "_", fleetnames[j]))
} else {
size_selex_label[[j]] <- paste0("SizeSel_", j, "PMalOff_", 1:nparms, "_", fleetnames[j])
}
size_selex_Nparms[j] <- size_selex_Nparms[j] + nparms
}
if (ctllist[["size_selex_types"]][j, 3] == 4) # has value 4 - differs by version
{
if (nver < 3.24) {
nparms <- 3
} else {
nparms <- 5
}
if (size_selex_Nparms[j] > 0) {
size_selex_label[[j]] <- c(size_selex_label[[j]], paste0("SizeSel_", j, "PFemOff_", 1:nparms, "_", fleetnames[j]))
} else {
size_selex_label[[j]] <- paste0("SizeSel_", j, "PFemOff_", 1:nparms, "_", fleetnames[j])
}
size_selex_Nparms[j] <- size_selex_Nparms[j] + nparms
}
}
age_selex_Nparms <- vector(mode = "numeric", length = Nfleet + Nsurveys)
age_selex_label <- list()
for (j in 1:(Nfleet + Nsurveys)) {
age_selex_Nparms[j] <- selex_patterns[age_selex_pattern_vec[j] + 1]
## some need special treatment
special_selex <- FALSE
if (age_selex_pattern_vec[j] == 17) { # random walk - check for maximum age in special
if (ctllist[["age_selex_types"]][["Special"]][j] > 0) {
special_selex <- TRUE
age_selex_Nparms[j] <- ctllist[["age_selex_types"]][["Special"]][j] + 1
age_selex_label[[j]] <- paste0("AgeSel_", j, "P_", 1:age_selex_Nparms[j], "_", fleetnames[j])
}
}
if (age_selex_pattern_vec[j] == 27) {
special_selex <- TRUE
age_selex_Nparms[j] <- age_selex_Nparms[j] + ctllist[["age_selex_types"]][j, 4] * 2
# age_selec_label[j]<-paste0("AgeSel_",j,"P_",1:agee_selex_Nparms[j],"_",fleetnames[j])
age_selex_label[[j]] <- c(
paste0("AgeeSpline_Code_", fleetnames[j], "_", j),
paste0("AgeSpline_GradLo_", fleetnames[j], "_", j),
paste0("AgeeSpline_GradHi_", fleetnames[j], "_", j),
paste0("AgeeSpline_Knot_", 1:ctllist[["age_selex_types"]][j, 4], "_", fleetnames[j], "_", j),
paste0("AgeeSpline_Val_", 1:ctllist[["age_selex_types"]][j, 4], "_", fleetnames[j], "_", j)
)
}
if (!special_selex) {
age_selex_label[[j]] <- if (age_selex_Nparms[j] > 0) {
paste0("AgeSel_", j, "P_", 1:age_selex_Nparms[j], "_", fleetnames[j])
} else {
NULL
}
}
# do extra retention parameters
if (ctllist[["age_selex_types"]][j, 2] > 0) # has discard type 1 or 2
{
if (age_selex_Nparms[j] > 0) {
age_selex_label[[j]] <- c(age_selex_label[[j]], paste0("AgeSel_", j, "PRet_", 1:4, "_", fleetnames[j]))
} else {
age_selex_label[[j]] <- paste0("AgeSel_", j, "PRet_", 1:4, "_", fleetnames[j])
}
age_selex_Nparms[j] <- age_selex_Nparms[j] + 4
}
if (ctllist[["age_selex_types"]][j, 2] == 2) # has discard type 2
{
if (age_selex_Nparms[j] > 0) {
age_selex_label[[j]] <- c(age_selex_label[[j]], paste0("AgeSel_", j, "PDis_", 1:4, "_", fleetnames[j]))
} else {
age_selex_label[[j]] <- paste0("AgeSel_", j, "PDis_", 1:4, "_", fleetnames[j])
}
age_selex_Nparms[j] <- age_selex_Nparms[j] + 4
}
# do extra offset parameters
if (ctllist[["age_selex_types"]][j, 3] == 1) # has value 1
{
if (age_selex_Nparms[j] > 0) {
age_selex_label[[j]] <- c(age_selex_label[[j]], paste0("AgeSel_", j, "PMale_", 1:4, "_", fleetnames[j]))
} else {
age_selex_label[[j]] <- paste0("AgeSel_", j, "PMale_", 1:4, "_", fleetnames[j])
}
age_selex_Nparms[j] <- age_selex_Nparms[j] + 4
}
if (ctllist[["age_selex_types"]][j, 3] == 2) # has value 2
{
if (age_selex_Nparms[j] > 0) {
age_selex_label[[j]] <- c(age_selex_label[[j]], paste0("AgeSel_", j, "PFemOff_", 1:4, "_", fleetnames[j]))
} else {
age_selex_label[[j]] <- paste0("AgeSel_", j, "PFemOff_", 1:4, "_", fleetnames[j])
}
age_selex_Nparms[j] <- age_selex_Nparms[j] + 4
}
if (ctllist[["age_selex_types"]][j, 3] == 3) # has value 3 - differs by version
{
if (nver < 3.24) {
nparms <- 3
} else {
nparms <- 5
}
if (age_selex_Nparms[j] > 0) {
age_selex_label[[j]] <- c(age_selex_label[[j]], paste0("AgeSel_", j, "PMalOff_", 1:nparms, "_", fleetnames[j]))
} else {
age_selex_label[[j]] <- paste0("AgeSel_", j, "PMalOff_", 1:nparms, "_", fleetnames[j])
}
age_selex_Nparms[j] <- age_selex_Nparms[j] + nparms
}
}
if (verbose) {
message("size_selex_Nparms\n")
print(size_selex_Nparms)
}
if (verbose) {
message("age_selex_Nparms\n")
print(age_selex_Nparms)
}
# Selex parlines ----
# _LO HI INIT PRIOR PR_type SD PHASE env-var use_dev dev_minyr dev_maxyr dev_stddev Block Block_Fxn
# Size selex
if (sum(size_selex_Nparms) > 0) {
ctllist <- add_df(ctllist,
name = "size_selex_parms", nrow = sum(size_selex_Nparms),
ncol = 14, col.names = lng_par_colnames,
comments = unlist(size_selex_label)
)
}
# Age selex
if (sum(age_selex_Nparms) > 0) {
ctllist <- add_df(ctllist,
name = "age_selex_parms", nrow = sum(age_selex_Nparms),
ncol = 14, col.names = lng_par_colnames,
comments = unlist(age_selex_label)
)
}
##########################
# selex timevarying ----
## Following pars are not yet implemented in this R code
# _Cond 0 #_custom_sel-env_setup (0/1)
# _Cond -2 2 0 0 -1 99 -2 #_placeholder when no enviro fxns
DoAdjust <- FALSE
if (sum(ctllist[["age_selex_parms"]][, 13]) + sum(ctllist[["size_selex_parms"]][, 13]) > 0) {
DoAdjust <- TRUE
ctllist <- add_elem(ctllist, "DoCustom_sel_blk_setup") # _custom_sel-blk_setup (0/1)
if (ctllist[["DoCustom_sel_blk_setup"]]) {
# FIND relevant blocks
blks <- ctllist[["blocks_per_pattern"]]
# SUBSET params for Block > 0
pbksa <- ctllist[["age_selex_parms"]][ctllist[["age_selex_parms"]][["Block"]] > 0, ]
pbkss <- ctllist[["size_selex_parms"]][ctllist[["size_selex_parms"]][["Block"]] > 0, ]
pbks <- rbind(pbksa, pbkss)
# FIND no. params per block
counts <- rle(sort(pbks[, "Block"]))
# MULTIPLY number of blocks by number of params per block
k0 <- sum(blks[counts[["values"]]] * counts[["lengths"]])
if (k0 > 0) {
ctllist <- add_df(ctllist,
name = "custom_sel_blk_setup", nrow = k0, ncol = 7,
col.names = srt_par_colnames
)
}
}
}
# _Cond No selex parm trends
if (sum(ctllist[["age_selex_parms"]][, 9]) + sum(ctllist[["size_selex_parms"]][, 9]) > 0) {
DoAdjust <- TRUE
ctllist <- add_elem(ctllist, name = "selparm_Dev_Phase") # selparm_Dev_Phase
} else {
# _Cond -4 # placeholder for selparm_Dev_Phase
}
ctllist[["DoAdjust"]] <- DoAdjust
if (DoAdjust) {
ctllist <- add_elem(ctllist, name = "selex_adjust_method")
# _env/block/dev_adjust_method (1=standard; 2=logistic trans to keep in base parm bounds; 3=standard w/ no bound check)
}
# _Cond 0 #_env/block/dev_adjust_method (1=standard; 2=logistic trans to keep in base parm bounds; 3=standard w/ no bound check)
#
# tagging ----
# Tag loss and Tag reporting parameters go next
ctllist <- add_elem(ctllist, name = "TG_custom") # TG_custom: 0=no read; 1=read if tags exist
if (ctllist[["TG_custom"]]) {
## Following parameters are to be read
## Fro each SS's tag release group
## . one initial Tag loss parameter
## . one continuos Tag loss parameter
## . NB over-dispersion paramater
## For each fleet
## . one tag reporting rate paramater
## . one tag reporting rate decay paramater
##
## In total N_tag_groups*3+ Nfleet*2 parameters are needed to read
##
# Initial tag loss
ctllist <- add_df(ctllist,
name = "TG_Loss_init", nrow = N_tag_groups, ncol = 14,
col.names = c(
srt_par_colnames,
"Dum", "Dum", "Dum", "Dum", "Dum", "Dum", "Dum"
)
)
# continuous tag loss
ctllist <- add_df(ctllist,
name = "TG_Loss_chronic", nrow = N_tag_groups, ncol = 14,
col.names = c(
srt_par_colnames,
"Dum", "Dum", "Dum", "Dum", "Dum", "Dum", "Dum"
),
comments = paste0("#_TG_Loss_chronic_", 1:N_tag_groups)
)
# NB over-dispersion
ctllist <- add_df(ctllist,
name = "TG_overdispersion", nrow = N_tag_groups, ncol = 14,
col.names = c(
srt_par_colnames,
"Dum", "Dum", "Dum", "Dum", "Dum", "Dum", "Dum"
),
comments = paste0("#_TG_overdispersion_", 1:N_tag_groups)
)
# TG_Report_fleet
ctllist <- add_df(ctllist,
name = "TG_Report_fleet", nrow = Nfleet, 14,
col.names = c(
srt_par_colnames,
"Dum", "Dum", "Dum", "Dum", "Dum", "Dum", "Dum"
)
)
# TG_Report_fleet_decay
ctllist <- add_df(ctllist,
name = "TG_Report_fleet_decay", nrow = Nfleet, 14,
col.names = c(
srt_par_colnames,
"Dum", "Dum", "Dum", "Dum", "Dum", "Dum", "Dum"
)
)
} else {
# _Cond -6 6 1 1 2 0.01 -4 0 0 0 0 0 0 0 #_placeholder if no parameters
#
}
# Var adj ----
ctllist <- add_elem(ctllist, "DoVar_adjust") # _Variance_adjustments_to_input_values
if (ctllist[["DoVar_adjust"]] > 0) {
ctllist <- add_df(ctllist,
name = "Variance_adjustments", nrow = 6, ncol = Nfleet + Nsurveys,
col.names = paste0("Fleet", 1:(Nfleet + Nsurveys))
)
rownames(ctllist[["Variance_adjustments"]]) <- paste0("#_", paste(c(
"add_to_survey_CV",
"add_to_discard_stddev",
"add_to_bodywt_CV",
"mult_by_lencomp_N",
"mult_by_agecomp_N",
"mult_by_size-at-age_N"
)))
}
# Lambdas ----
ctllist <- add_elem(ctllist, "maxlambdaphase") # _maxlambdaphase
ctllist <- add_elem(ctllist, "sd_offset") # _sd_offset
ctllist <- add_elem(ctllist, "N_lambdas") # number of changes to make to default Lambdas (default value is 1.0)
if (ctllist[["N_lambdas"]] > 0) {
ctllist <- add_df(ctllist,
name = "lambdas", nrow = ctllist[["N_lambdas"]], ncol = 5,
col.names = c("like_comp", "fleet/survey", "phase", "value", "sizefreq_method")
)
# find and delete duplicates
chk1 <- duplicated(ctllist[["lambdas"]])
if (any(chk1)) # there are duplicates
{
ctllist[["lambdas"]] <- ctllist[["lambdas"]][!chk1, ]
ctllist[["N_lambdas"]] <- nrow(ctllist[["lambdas"]])
ctllist[["warnings"]] <- paste(ctllist[["warnings"]], "Duplicate_lambdas", sep = ",")
}
for (i in 1:ctllist[["N_lambdas"]]) {
like_comp <- ctllist[["lambdas"]][i, 1]
fl <- ctllist[["lambdas"]][i, 2]
phz <- ctllist[["lambdas"]][i, 3]
value <- ctllist[["lambdas"]][i, 4]
sizefreq_method <- ctllist[["lambdas"]][i, 5]
if (like_comp == 1) {
rownames(ctllist[["lambdas"]])[i] <- paste0("Surv_", fleetnames[fl], "_Phz", phz)
} else if (like_comp == 4) {
rownames(ctllist[["lambdas"]])[i] <- paste0("length_", fleetnames[fl], "_sizefreq_method_", sizefreq_method, "_Phz", phz)
} else if (like_comp == 5) {
rownames(ctllist[["lambdas"]])[i] <- paste0("age_", fleetnames[fl], "_Phz", phz)
} else if (like_comp == 5) {
sizefreq_method <- ctllist[["lambdas"]][i, 5]
rownames(ctllist[["lambdas"]])[i] <- paste0("SizeFreq_", fleetnames[fl], "_sizefreq_method_", sizefreq_method, "_Phz", phz)
} else if (like_comp == 8) {
rownames(ctllist[["lambdas"]])[i] <- paste0("catch_", fleetnames[fl], "_Phz", phz)
} else if (like_comp == 9) {
rownames(ctllist[["lambdas"]])[i] <- paste0("init_equ_catch_", fleetnames[fl], "_lambda_for_init_equ_catch_can_only_enable/disable for_all_fleets_Phz", phz)
} else if (like_comp == 10) {
rownames(ctllist[["lambdas"]])[i] <- paste0("recrdev_Phz", phz)
} else if (like_comp == 11) {
rownames(ctllist[["lambdas"]])[i] <- paste0("parm_prior_", fleetnames[fl], "_Phz", phz)
} else if (like_comp == 12) {
rownames(ctllist[["lambdas"]])[i] <- paste0("parm_dev_Phz", phz)
} else if (like_comp == 13) {
rownames(ctllist[["lambdas"]])[i] <- paste0("CrashPen_Phz", phz)
} else if (like_comp == 15) {
rownames(ctllist[["lambdas"]])[i] <- paste0("Tag-comp-likelihood-", fl, "_Phz", phz)
} else if (like_comp == 16) {
rownames(ctllist[["lambdas"]])[i] <- paste0("Tag-negbin-likelihood-", fl, "_Phz", phz)
} else if (like_comp == 17) {
rownames(ctllist[["lambdas"]])[i] <- paste0("F-ballpark-", fl, "_Phz", phz)
}
}
}
# Like_comp codes: 1=surv; 2=disc; 3=mnwt; 4=length; 5=age; 6=SizeFreq; 7=sizeage; 8=catch;
# 9=init_equ_catch; 10=recrdev; 11=parm_prior; 12=parm_dev; 13=CrashPen; 14=Morphcomp; 15=Tag-comp; 16=Tag-negbin
# more sd reporting ----
ctllist <- add_elem(ctllist, "more_stddev_reporting") # (0/1) read specs for more stddev reporting
if (ctllist[["more_stddev_reporting"]] != 0) {
ctllist <- add_vec(ctllist, name = "stddev_reporting_specs", length = 9)
## Selex bin
if (ctllist[["stddev_reporting_specs"]][4] > 0) {
ctllist <-
add_vec(ctllist, name = "stddev_reporting_selex", length = ctllist[["stddev_reporting_specs"]][4])
}
## Growth bin
if (ctllist[["stddev_reporting_specs"]][6] > 0) {
ctllist <-
add_vec(ctllist, name = "stddev_reporting_growth", length = ctllist[["stddev_reporting_specs"]][6])
}
## N at age
if (ctllist[["stddev_reporting_specs"]][9] > 0) {
ctllist <-
add_vec(ctllist, name = "stddev_reporting_N_at_A", length = ctllist[["stddev_reporting_specs"]][9])
}
}
if (ctllist$".dat"[ctllist$".i"] == 999) {
if (verbose) message("read of control file complete (final value = 999)")
ctllist[["eof"]] <- TRUE
} else {
message(
"Error: final value is", ctllist$".dat"[ctllist$".i"],
", but should be 999\n"
)
ctllist[["eof"]] <- FALSE
}
ctllist$".dat" <- NULL
ctllist$".i" <- NULL
# return the result
return(ctllist)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.