mw_ss.R

library(r4ss)  # NB: This is to force shinyapps to register the dependency
# https://r4ss.github.io/r4ss/articles/r4ss-intro-vignette.html
# https://nmfs-stock-synthesis.github.io/doc/SS330_User_Manual_release.html

template_str <- function (s) {
    stringr::str_interp(s, parent.frame(1))
}
escape_sym <- Vectorize(function (s) deparse1(as.symbol(s), backtick = TRUE))

mw_ss_code_readxl <- function (sheet_name, xlsx) {
    if (!nzchar(xlsx)) return("")
    template_str('${escape_sym(sheet_name)} <- readxl::read_excel(data_path, ${deparse1(sheet_name)}, na = c("", "NA"))\n')
}

mw_ss_code_header <- function (spec, xlsx) {
    mod_dir <- gsub('\\.xlsx$' , '-ss3', xlsx)

    template_str(r'(
library(r4ss)
library(reshape2)
library(readxl)
library(dplyr)

mod_path <- ${deparse1(mod_dir)}
data_path <- ${deparse1(xlsx)}

# parse_levels from gadget3
parse_levels <- function (lvls, var_name) {
    m <- suppressWarnings(as.numeric(lvls))
    if (!anyNA(m)) return(data.frame(
        names = lvls,
        lower_incl = TRUE,
        lower_bound = m,
        upper_bound = c(tail(m, -1), Inf),  # NB: No data about final bound, assume open-ended
        upper_incl = FALSE,
        open_ended_upper = TRUE,
        stringsAsFactors = FALSE))

    m <- regmatches(lvls, regexec("^(\\[|\\()(.*),(.*)(\\]|\\))", lvls))
    if (all(vapply(m, length, numeric(1)) == 5)) return(data.frame(
        names = lvls,
        lower_incl = vapply(m, function (mm) identical(mm[[2]], '['), logical(1)),
        lower_bound = vapply(m, function (mm) as.numeric(mm[[3]]), numeric(1)),
        upper_bound = vapply(m, function (mm) as.numeric(mm[[4]]), numeric(1)),
        upper_incl = vapply(m, function (mm) identical(mm[[5]], ']'), logical(1)),
        open_ended_upper = is.infinite(as.numeric(tail(m, 1)[[1]][[4]])),
        stringsAsFactors = FALSE))

    stop("Unknown form of ", var_name, " levels, see ?cut for formatting: ", paste(lvls, collapse = ", "))
}

# Create new SS3 model based on "simple_small"
r4ss::copy_SS_inputs(
    dir.old = system.file(file.path("extdata", "simple_small"), package = "r4ss"),
    dir.new = mod_path,
    overwrite = TRUE)
inputs <- r4ss::SS_read(dir = mod_path)

inputs$dat$Comments <- sprintf("#C %s model generated by ModelWizard", mod_path)

# Clear tables ####################
## data file ----
inputs$dat$fleetinfo <- inputs$dat$fleetinfo[c(),,drop = FALSE]
inputs$dat$catch     <- inputs$dat$catch[c(),,drop = FALSE]
inputs$dat$CPUEinfo  <- inputs$dat$CPUEinfo[c(),,drop = FALSE]
inputs$dat$CPUE      <- inputs$dat$CPUE[c(),,drop = FALSE]
inputs$dat$len_info  <- inputs$dat$len_info[c(),,drop = FALSE]
inputs$dat$lencomp   <- inputs$dat$lencomp[c(),,drop = FALSE]
inputs$dat$age_info  <- inputs$dat$age_info[c(),,drop = FALSE]
inputs$dat$agecomp   <- inputs$dat$agecomp[c(),,drop = FALSE]
inputs$dat$ageerror  <- inputs$dat$ageerror[c(),,drop = FALSE]
inputs$dat$MeanSize_at_Age_obs <- inputs$dat$MeanSize_at_Age_obs[c(),,drop = FALSE]
## control file ----
inputs$ctl$MG_parms         <- inputs$ctl$MG_parms[c(),,drop = FALSE]
inputs$ctl$SR_parms         <- inputs$ctl$SR_parms[c(),,drop = FALSE]
inputs$ctl$Q_options        <- inputs$ctl$Q_options[c(),,drop = FALSE]
inputs$ctl$Q_parms          <- inputs$ctl$Q_parms[c(),,drop = FALSE]
inputs$ctl$size_selex_types <- inputs$ctl$size_selex_types[c(),,drop = FALSE]
inputs$ctl$size_selex_parms <- inputs$ctl$size_selex_parms[c(),,drop = FALSE]
inputs$ctl$age_selex_types  <- inputs$ctl$age_selex_types[c(),,drop = FALSE]
inputs$ctl$age_selex_parms  <- inputs$ctl$age_selex_parms[c(),,drop = FALSE]
inputs$ctl$lambdas          <- inputs$ctl$lambdas[c(),,drop = FALSE]
)')}

mw_ss_code_area <- function (spec) {
    area_count <- length(spec$area)
    template_str(r'(
# Create area definitions ####################
inputs$dat$N_areas <- ${deparse1(area_count)}
)')}

mw_ss_code_time <- function (r, spec) {
    year_max <- as.integer(r$year_max)
    year_min <- as.integer(r$year_min)
    steps_count <- as.integer(r$steps)
    step_lengths <- rep(as.integer(12 / steps_count), steps_count)

    template_str(r'(
# Create time definitions ####################

inputs$dat$styr <- ${deparse1(year_min)}

inputs$dat$endyr <- ${deparse1(year_max)}
inputs$dat$nseas <- ${deparse1(steps_count)}
inputs$dat$months_per_seas <- rep(12/inputs$dat$nseas, inputs$dat$nseas)
)')}

mw_ss_code_stock <- function (r, spec, xlsx) {
    if (nrow(spec$stock) > 1) stop("Only one stock supported, got: ", paste(spec$stock$name))

    template_str(r'(
# Create stock definition for ${r$name} ####################

inputs$dat$spawn_month <- 1
inputs$dat$Nsexes <- 1
inputs$dat$Nages <- ${deparse1(r$age_max - r$age_min)}
inputs$dat$lbin_method  <-2 # option SS3
inputs$dat$binwidth     <- ${deparse1(r$lg_size)}
inputs$dat$minimum_size <- ${deparse1(r$lg_min)}
inputs$dat$maximum_size <- ${deparse1(r$lg_max)}
inputs$dat$use_lencomp  <-1 # option SS3
inputs$dat$lbin_vector <- seq(${deparse1(r$lg_min)}, ${deparse1(r$lg_max)}, ${deparse(r$lg_size)})
inputs$dat$agebin_vector <- seq(${deparse1(r$age_min)}, ${deparse1(r$age_max)}, 1)

# Growth
inputs$ctl$Growth_Age_for_L1<-0
inputs$ctl$Growth_Age_for_L2<-999
inputs$ctl$MG_parms["L_at_Amin", c("LO","HI","INIT","PHASE")]  <- c(0,10, 5,-3)
inputs$ctl$MG_parms["L_at_Amax", c("LO","HI","INIT","PHASE")]  <- c(5,22,19,-3)
inputs$ctl$MG_parms["VonBert_k", c("LO","HI","INIT","PHASE")]  <- c(0.1,2,0.89,-3)
inputs$ctl$MG_parms["CV_young", c("LO","HI","INIT","PHASE")]   <- c(0.05,0.25,0.1,-3)
inputs$ctl$MG_parms["CV_old", c("LO","HI","INIT","PHASE")]     <- c(0.05,0.25,0.1,-3)
inputs$ctl$MG_parms["Wtlen_1", c("LO","HI","INIT","PHASE")]    <- c(-3,3,0.00000313,-3) # kg
inputs$ctl$MG_parms["Wtlen_2", c("LO","HI","INIT","PHASE")]    <- c(-3,3,3.278,-3) # kg
inputs$ctl$MG_parms["Mat50%", c("LO","HI","INIT","PHASE")]     <- c(-3,15,22,-3)
inputs$ctl$MG_parms["Mat_slope", c("LO","HI","INIT","PHASE")]  <- c(-3,3,-0.45,-3)
inputs$ctl$MG_parms["Eggs_alpha", c("LO","HI","INIT","PHASE")] <- c(-3,3,1,-3)
inputs$ctl$MG_parms["Eggs_beta", c("LO","HI","INIT","PHASE")]  <- c(-3,3,0,-3)
inputs$ctl$MG_parms["CohortGrowDev", c("LO","HI","INIT","PHASE")] <- c(1,1,1,-3)
inputs$ctl$MG_parms["FracFemale", c("LO","HI","INIT","PHASE")]  <- c(0.0000001,0.9999999,0.5,-3)
inputs$ctl$MG_parms[is.na(inputs$ctl$MG_parms)]<-0

# Stock-recruitment (SR)
inputs$ctl$SR_function <- 4
inputs$ctl$SR_parms["SR_LN(R0)", c("LO","HI","INIT","PHASE")] <- c(5,10,9,1)
inputs$ctl$SR_parms["SR_SCAA_null", c("LO","HI","INIT","PHASE")] <- c(0.2,1,0.9,-4)
inputs$ctl$SR_parms["SR_sigmaR", c("LO","HI","INIT","PHASE")] <- c(0,2,0.6,-4)
inputs$ctl$SR_parms["SR_regime", c("LO","HI","INIT","PHASE")] <- c(-5,5,0,-4)
inputs$ctl$SR_parms["SR_autocorr", c("LO","HI","INIT","PHASE")] <- c(0,0,0,-4)
inputs$ctl$SR_parms[is.na(inputs$ctl$SR_parms)]<-0

# Natural Mortality
inputs$ctl$natM_type <- 3  # Age-varying
inputs$ctl$natM <- as.data.frame(sapply(
    paste("Age", inputs$dat$agebin_vector, sep = "_"),
    function (x) list(1.3) ))
rownames(inputs$ctl$natM) <- paste0("natM", seq_len(nrow(inputs$ctl$natM)))

# Recruitment deviation (recDev)
inputs$ctl$MainRdevYrFirst<-inputs$dat$styr
inputs$ctl$MainRdevYrLast<-inputs$dat$endyr
inputs$ctl$recdev_phase<-1
inputs$ctl$recdev_adv<-0
)')}


mw_ss_code_comm <- function (r, spec, xlsx) {
    fleet_sym <- escape_sym(paste0('fleet_', r$name))
    area_names <- spec$area$name
    stock_list <- lapply(spec$stock$name, as.symbol)
    data_name <- paste("landings", r$name, sep = "_")
    data_sym <- escape_sym(data_name)

    if (r$step_active == 0) {
        surveytiming <- "-1" # All steps
    } else {
        surveytiming <- paste0("(", r$step_active, "* 12/inputs$dat$nseas)/12")
    }

    template_str(r'(
# Create commercial fleet ${r$name} ####################
inputs$dat$fleetinfo <- rbind(inputs$dat$fleetinfo, data.frame(
    type = 1,
    surveytiming = ${surveytiming},
    area = 1,  # Assume one area
    units = ${if (r$landings == "weight") "1,  # Biomass (metric tons)" else "2,  # Numbers (thousands of fish)"}
    need_catch_mult = 0,
    fleetname = ${deparse1(r$name)},
    stringsAsFactors = FALSE))

inputs$dat$len_info[${deparse1(r$name)},] <- list(
    mintailcomp = -1,
    addtocomp = 1e-4,
    combine_M_F = 0,
    CompressBins = 0,
    CompError = 0,
    ParmSelect = 0,
    minsamplesize = 1)

inputs$dat$age_info[${deparse1(r$name)},] <- list(
    mintailcomp = -1,
    addtocomp = 1e-4,
    combine_M_F = 0,
    CompressBins = 0,
    CompError = 0,
    ParmSelect = 0,
    minsamplesize = 1)

inputs$dat$CPUEinfo[${deparse1(r$name)},] <- list(
    Fleet = which(inputs$dat$fleetinfo$fleetname == ${deparse1(r$name)}),
    Units = 1,
    Errtype = 0,
    SD_Report = 0)

# Equilibrium catch fields
if (nrow(inputs$dat$catch) == 0) inputs$dat$catch <- rbind(inputs$dat$catch, data.frame(
    year = -999,
    seas = seq_len(inputs$dat$nseas),
    fleet = which(inputs$dat$fleetinfo$fleetname == ${deparse1(r$name)}),
    catch = 0,
    catch_se = 0.1,
    stringsAsFactors = TRUE))

${mw_ss_code_readxl(data_name, xlsx)}
inputs$dat$catch <- rbind(inputs$dat$catch, data.frame(
    year = ${data_sym}$year,
    seas = ${data_sym}$step,
    fleet = which(inputs$dat$fleetinfo$fleetname == ${deparse1(r$name)}),
    catch = ${data_sym}[[${deparse1(r$landings)}]],
    catch_se = 0.1,
    stringsAsFactors = TRUE))

inputs$ctl$size_selex_types[${deparse1(r$name)}, "Pattern"] <- 1
inputs$ctl$size_selex_types[is.na(inputs$ctl$size_selex_types)] <- 0

inputs$ctl$size_selex_parms[paste0("SizeSel_P_1_", ${deparse1(r$name)}),c("LO","HI","INIT","PHASE")] <- c(-1,10,8,2)
inputs$ctl$size_selex_parms[paste0("SizeSel_P_2_", ${deparse1(r$name)}),c("LO","HI","INIT","PHASE")] <- c(-1,20,10,2)
inputs$ctl$size_selex_parms[is.na(inputs$ctl$size_selex_parms)] <- 0

inputs$ctl$age_selex_types[${deparse1(r$name)}, "Pattern"] <- 12
inputs$ctl$age_selex_types[is.na(inputs$ctl$age_selex_types)] <- 0

inputs$ctl$age_selex_parms[paste0("AgeSel_P_1_", ${deparse1(r$name)}),c("LO","HI","INIT","PHASE")] <- c(-2,5.5,0.1,-1)
inputs$ctl$age_selex_parms[paste0("AgeSel_P_2_", ${deparse1(r$name)}),c("LO","HI","INIT","PHASE")] <- c(-1,5.5,5,-1)
inputs$ctl$age_selex_parms[is.na(inputs$ctl$age_selex_parms)] <- 0
)')}

mw_ss_code_surv <- function (r, spec, xlsx) {
    fleet_sym <- escape_sym(r$name)
    stock_list <- lapply(spec$stock$name, as.symbol)

    if (r[['si']] == 'none') return("")

    cpue_name <- unname(paste('si', r$name, sep = "_"))
    cpue_sym <- escape_sym(cpue_name)

    if (r$step_active == 0) {
        surveytiming <- "-1" # All steps
    } else {
        surveytiming <- paste0("(", r$step_active, "* 12/inputs$dat$nseas)/12")
    }

    template_str(r'(
# Create survey fleet ${r$name} ####################
${mw_ss_code_readxl(cpue_name, xlsx)}

inputs$dat$fleetinfo <- rbind(inputs$dat$fleetinfo, data.frame(
    type = 3,
    surveytiming = ${surveytiming},
    area = 1,  # Assume one area
    units = 1,  # NB: Survey fleets ignore this value
    need_catch_mult = 0,
    fleetname = ${deparse1(r$name)},
    stringsAsFactors = FALSE))

inputs$dat$len_info[${deparse1(r$name)},] <- list(
    mintailcomp = -1,
    addtocomp = 1e-4,
    combine_M_F = 0,
    CompressBins = 0,
    CompError = 0,
    ParmSelect = 0,
    minsamplesize = 1)

inputs$dat$age_info[${deparse1(r$name)},] <- list(
    mintailcomp = -1,
    addtocomp = 1e-4,
    combine_M_F = 0,
    CompressBins = 0,
    CompError = 0,
    ParmSelect = 0,
    minsamplesize = 1)

inputs$dat$CPUEinfo[${deparse1(r$name)},] <- list(
    Fleet = which(inputs$dat$fleetinfo$fleetname == ${deparse1(r$name)}),
    Units = ${if (r$si == "weight") "1,  # Biomass (metric tons)" else "2,  # Numbers (thousands of fish)"},
    Errtype = 0,
    SD_Report = 0)

inputs$dat$CPUE <- rbind(inputs$dat$CPUE, data.frame(
  year = ${cpue_sym}$year,
  seas = ${cpue_sym}$step*(12/inputs$dat$nseas),
  fleet = which(inputs$dat$fleetinfo$fleetname == ${deparse1(r$name)}),
  catch = ${cpue_sym}[["weight"]],
  catch_se = 0.3,
  stringsAsFactors = TRUE))

# https://nmfs-ost.github.io/ss3-doc/SS330_User_Manual_release.html#catchability CPUE scaling
inputs$ctl$Q_options[${deparse1(r$name)},] <- list(
    fleet = which(inputs$dat$fleetinfo$fleetname == ${deparse1(r$name)}),
    link = 1,
    link_info = 0,
    extra_se = 0,
    biasadj = 0,
    float = 0 )

inputs$ctl$Q_parms[paste0("LnQ_base_", ${deparse1(r$name)}),] <- list(
    LO = 0,
    HI = 5.5,
    INIT = 0.1,
    PRIOR = 0,
    PR_SD = 0,
    PR_type = 0,
    PHASE = 2,
    "env-var" = 0,
    dev_mnyr = 0,
    dev_mxyr = 0,
    dev_PH = 0,
    Block = 0,
    Blk_Fxn = 0)

inputs$ctl$size_selex_types[${deparse1(r$name)}, "Pattern"] <- 1
inputs$ctl$size_selex_types[is.na(inputs$ctl$size_selex_types)] <- 0

inputs$ctl$size_selex_parms[paste0("SizeSel_P_1_", ${deparse1(r$name)}),c("LO","HI","INIT","PHASE")] <- c(-1,15,10,3)
inputs$ctl$size_selex_parms[paste0("SizeSel_P_2_", ${deparse1(r$name)}),c("LO","HI","INIT","PHASE")] <- c(-1,20,12,3)
inputs$ctl$size_selex_parms[is.na(inputs$ctl$size_selex_parms)] <- 0

inputs$ctl$age_selex_types[${deparse1(r$name)}, "Pattern"] <- 12
inputs$ctl$age_selex_types[is.na(inputs$ctl$age_selex_types)] <- 0

inputs$ctl$age_selex_parms[paste0("AgeSel_P_1_", ${deparse1(r$name)}),c("LO","HI","INIT","PHASE")] <- c(-2,5.5,0.1,-1)
inputs$ctl$age_selex_parms[paste0("AgeSel_P_2_", ${deparse1(r$name)}),c("LO","HI","INIT","PHASE")] <- c(-1,5.5,5,-1)
inputs$ctl$age_selex_parms[is.na(inputs$ctl$age_selex_parms)] <- 0
)')}

mw_ss_code_ldist <- function (r, spec, xlsx) {
    fleet_sym <- escape_sym(r$name)
    stock_list <- lapply(spec$stock$name, as.symbol)

    if (r[['ldist']] == 'none') return("")

    ldist_name <- unname(paste('ldist', r$name, sep = "_"))
    ldist_sym <- escape_sym(ldist_name)
    
    template_str(r'(
# Create length distribution for ${r$name} ####################

${mw_ss_code_readxl(ldist_name, xlsx)}

if (nrow(inputs$dat$lencomp) == 0) {
    # Seed lencomp with columns based on lbin_vector
    cols <- paste0("L", inputs$dat$lbin_vector)
    cols <- structure(as.list(rep(0, length(cols))), names = cols)
    inputs$dat$lencomp <- as.data.frame(c(list(
        Yr = 1999,
        Seas = 1,
        Fleet = 0,
        Gender = 0,
        Part = 0,
        NSamp = 100), cols))[c(),]
}

# Convert lengths into SS-compatible names
${ldist_sym}$length <- as.factor(${ldist_sym}$length)
lvls <- parse_levels(levels(${ldist_sym}$length))
lvls_map <- structure(paste0("L", lvls$lower_bound), names = lvls$names)
levels(${ldist_sym}$length) <- lvls_map[levels(${ldist_sym}$length)]

# Rotate length groupings, renaming to SS labels
${ldist_sym} <- reshape2::dcast(
    ${ldist_sym},
    year + step + area ~ length,
    value.var = ${deparse1(r$ldist)})

inputs$dat$lencomp <- dplyr::bind_rows(inputs$dat$lencomp, cbind(data.frame(
    Yr = ${ldist_sym}$year,
    Seas = ${ldist_sym}$step * (12/inputs$dat$nseas),
    Fleet = which(inputs$dat$fleetinfo$fleetname == ${deparse1(r$name)}),
    Gender = 0,
    Part = 0,
    NSamp = 100,
    stringsAsFactors = TRUE), ${ldist_sym}[,lvls_map]))
)')}

mw_ss_code_aldist <- function (r, spec, xlsx) {
    fleet_sym <- escape_sym(r$name)
    stock_list <- lapply(spec$stock$name, as.symbol)

    if (r[['aldist']] == 'none') return("")

    data_name <- paste("aldist", r$name, sep = "_")
    data_sym <- escape_sym(data_name)

    agecomp_sym <- escape_sym(paste("agecomp", r$name, sep = "_"))
    agedist_sym <- escape_sym(paste("agedist", r$name, sep = "_"))
    agelencomp_sym <- escape_sym(paste("agelencomp", r$name, sep = "_"))

    template_str(r'(
# Create age-length distribution for ${r$name} ####################

if (nrow(inputs$dat$agecomp) == 0) {
  # Seed agecomp with columns based on abin_vector
  cols <- paste0("E", inputs$dat$agebin_vector)
  cols <- structure(as.list(rep(0, length(cols))), names = cols)

  inputs$dat$agecomp <- as.data.frame(c(list(
    Yr = 1999,
    Seas = 1,
    Fleet = 1,
    Gender = 0,
    Part = 0,
    Ageerr= 1,
    Lbin_lo=0,
    Lbin_hi=0,
    NSamp = 100), cols))[c(),]
}

${mw_ss_code_readxl(data_name, xlsx)}

${data_sym}$age <- as.factor(${data_sym}$age)
lvls_a <- parse_levels(levels(${data_sym}$age))
lvls_a_map <- structure(paste0("E", lvls_a$lower_bound), names = lvls_a$names)
levels(${data_sym}$age) <- lvls_a_map[levels(${data_sym}$age)]

inputs$dat$ageerror <- as.data.frame(t(data.frame(
    "1" = inputs$dat$agebin_vector + .5,
    "2" = 0.001,
    row.names = paste0("age", inputs$dat$agebin_vector))))

# *Age composition*
${agedist_sym}<-${data_sym} |>
  group_by(year,step,age) |> summarise(total=sum(number, na.rm = TRUE))
# Rotate age groupings, renaming to SS labels
${agedist_sym} <- reshape2::dcast(${agedist_sym},
                                year + step  ~ age, value.var = "total")

${agecomp_sym}<-cbind(data.frame(
  Yr = ${agedist_sym}$year,
  Seas = ${agedist_sym}$step*(12/inputs$dat$nseas),
  Fleet = -which(inputs$dat$fleetinfo$fleetname == ${deparse1(r$name)}),
  Gender = 0,
  Part = 0,
  Ageerr= 1,
  Lbin_lo=-1,
  Lbin_hi=-1,
  NSamp = rowSums(${agedist_sym}[,lvls_a_map], na.rm = TRUE),
  stringsAsFactors = TRUE), ${agedist_sym}[,lvls_a_map])
${agecomp_sym}[,lvls_a_map]<-${agecomp_sym}[,lvls_a_map]/${agecomp_sym}$NSamp #replace number by proportion

# *Age at length*
# Rotate age at length groupings, renaming to SS labels
${data_sym} <- reshape2::dcast(${data_sym},
  year + step + area + length ~ age, value.var = "number")

${data_sym}$length<-as.factor(${data_sym}$length)
lvls_l <- parse_levels(levels(${data_sym}$length))
lvls_l_map <- structure(lvls_l$lower_bound, names = lvls_l$names)
levels(${data_sym}$length) <- lvls_l_map[levels(${data_sym}$length)]

${agelencomp_sym}<-cbind(data.frame(
  Yr = ${data_sym}$year,
  Seas = ${data_sym}$step*(12/inputs$dat$nseas),
  Fleet = which(inputs$dat$fleetinfo$fleetname == ${deparse1(r$name)}),
  Gender = 0,
  Part = 0,
  Ageerr= 1,
  Lbin_lo=as.numeric(${data_sym}$length),
  Lbin_hi=as.numeric(${data_sym}$length),
  NSamp = rowSums(${data_sym}[,lvls_a_map], na.rm = TRUE),
  stringsAsFactors = TRUE), ${data_sym}[,lvls_a_map])
${agelencomp_sym}<-${agelencomp_sym}[!${agelencomp_sym}$NSamp==0,] # remove row with  size sample equal to 0.
${agelencomp_sym}[,lvls_a_map]<-${agelencomp_sym}[,lvls_a_map]/${agelencomp_sym}$NSamp #replace number by proportion

inputs$dat$agecomp <- dplyr::bind_rows(inputs$dat$agecomp, ${agecomp_sym}, ${agelencomp_sym})
inputs$dat$agecomp[is.na(inputs$dat$agecomp)] <- 0
rownames(inputs$dat$agecomp) <- NULL
)')}

mw_ss_code_params <- function (r, spec, xlsx) {
    stock_name <- spec$stock$name

    template_str(r'(
# Set params ###############################

inputs$ctl$MG_parms[grepl("VonBert_K", rownames(inputs$ctl$MG_parms)), "INIT"] <- 0.3
inputs$ctl$MG_parms[grepl("VonBert_K", rownames(inputs$ctl$MG_parms)), "LO"] <- 0.04
inputs$ctl$MG_parms[grepl("VonBert_K", rownames(inputs$ctl$MG_parms)), "HI"] <- 1.2

## Fishing Mortality (F) ----

# Setting the Ballpark Year for Fishing Mortality:
#   F_ballpark_year is set to a negative value of the end year (endyr) from the input data.
#   This indicates that the fishing mortality parameter is estimated or initialized relative to
#   the end year of the data.
# Initializing Fishing Mortality:
#   init_F is set to NULL, indicating that the initial fishing mortality values are not
#.  explicitly specified in this script. They might be initialized or estimated through other means in the model.
inputs$ctl$F_ballpark_year<- -inputs$dat$endyr
inputs$ctl$init_F<-NULL

)')}

mw_ss_code_footer <- function(spec, xlsx) {
    template_str(r'(
## Lambdas ----

# Lambdas:
#   maxlambdaphase is set to 1, indicating that lambda values will be estimated in
#.  the first phase of the model.
# The number of lambdas (N_lambdas) is determined based on the number of rows in lambdas.
#.If lambdas is empty, itss set to NULL.
# more_stddev_reporting is set to 0, indicating that no further standard deviation reporting
# is required.

inputs$ctl$maxlambdaphase<- 1

inputs$ctl$N_lambdas <- nrow(inputs$ctl$lambdas)
inputs$ctl$more_stddev_reporting<-0

# Finalise model #######################

# Set rownames (so they can be added as comments into data.ss)
rownames(inputs$dat$fleetinfo) <- inputs$dat$fleetinfo$fleetname

# Zero NA values in catch
for (n in names(inputs$dat$lencomp)[-seq_len(6)]) {
    inputs$dat$lencomp[[n]][is.na(inputs$dat$lencomp[[n]])] <- 0
}

# Derived values
inputs$dat$Nfleets <- sum(inputs$dat$fleetinfo$type %in% c(1,3))
inputs$dat$N_lbins <- length(inputs$dat$lbin_vector)
if (nrow(inputs$dat$lencomp) > 0) {
    inputs$dat$use_lencomp <- 1
} else {
    inputs$dat$use_lencomp <- 0
    inputs$dat$lencomp <- NULL  # r4ss wont write 0-row data.frames
}
inputs$dat$N_ageerror_definitions <- floor(nrow(inputs$dat$ageerror) / 2)
inputs$dat$N_agebins <- length(inputs$dat$agebin_vector)
if (nrow(inputs$dat$MeanSize_at_Age_obs) > 0) {
    inputs$dat$use_MeanSize_at_Age_obs <- 1
} else {
    inputs$dat$use_MeanSize_at_Age_obs <- 0
    inputs$dat$MeanSize_at_Age_obs <- NULL  # r4ss wont write 0-row data.frames
}
inputs$dat$Ngenders <- inputs$dat$Nsexes

inputs$ctl$fleetnames <- inputs$dat$fleetinfo$fleetname
inputs$ctl$Nfleets <- length(inputs$ctl$fleetnames)
inputs$ctl$N_lambdas <- nrow(inputs$ctl$lambdas)
if (nrow(inputs$ctl$lambdas) == 0) inputs$ctl$lambdas <- NULL
inputs$ctl$maxlambdaphase<- 1

# *------------------------------------------------------*
#  RUN AND OUTPUTS

# Write back to mod_path ----
# 1. Write Back to mod_path:
#    The SS_write function from the r4ss package is used to write the inputs back to the model path
#    specified by mod_path. If files with the same names already exist in the directory,
#    they will be overwritten (overwrite = TRUE).

r4ss::SS_write(inputs, dir = mod_path, overwrite = TRUE)

# Run model ----
# 2. Run Model:
#   The script changes the directory to the model path using the system function in the terminal.
#   File permissions are set to allow the computer to execute the models executable file (ss3).
#   The model SS3 is executed using the run function from the r4ss package, specifying the model directory,
#   executable file name (ss3), and other options like skipfinished and show_in_console.

ss3_exe <- r4ss::get_ss3_exe()
r4ss::run(dir = mod_path, exe = ss3_exe, skipfinished = FALSE, show_in_console = TRUE)

# Read the model output and print diagnostic messages ----
# 3. Read Model Output and Print Diagnostic Messages:
#    The SS_output function from the r4ss package is used to read the model output from the directory
#    specified by mod_path.
#    Verbose mode is turned on (verbose = TRUE) to print diagnostic messages.
#    Summary statistics are printed (printstats = TRUE) to provide insights into the models performance.

replist <- SS_output(
  dir = mod_path,
  verbose = TRUE,
  printstats = TRUE
)

# plots the results ----
# 4. Plot Results:
#    Finally, the SS_plots function from the r4ss package is used to generate plots based on the
#.   model output (replist).
  SS_plots(replist)
)')}

mw_ss_script <- function (
        spec,
        xlsx = "",
        compile = FALSE,
        run = FALSE) {
    stopifnot(is.list(spec) || is.environment(spec))
    stopifnot(length(intersect(names(spec), c("area", "comm", "surv", "stock", "time"))) == 5)

    # Run fn(row, ...) for each row in tbl
    row_apply <- function (tbl, fn, ...) vapply(
        seq_len(nrow(tbl)),
        function (i) fn(as.list(tbl[i,, drop = FALSE]), ...),
        character(1))

    paste(c(
        mw_ss_code_header(spec, xlsx),
        mw_ss_code_area(spec),
        row_apply(spec$time, mw_ss_code_time, spec),
        row_apply(spec$stock, mw_ss_code_stock, spec, xlsx),
        row_apply(spec$comm, mw_ss_code_comm, spec, xlsx),
        row_apply(spec$surv, mw_ss_code_surv, spec, xlsx),
        row_apply(spec$comm, mw_ss_code_ldist, spec, xlsx),
        row_apply(spec$surv, mw_ss_code_ldist, spec, xlsx),
        row_apply(spec$comm, mw_ss_code_aldist, spec, xlsx),
        row_apply(spec$surv, mw_ss_code_aldist, spec, xlsx),
        mw_ss_code_params(spec$params, spec, xlsx),
        mw_ss_code_footer(spec, xlsx),
#        (if (compile) mw_ss_code_compile(spec, xlsx) else ""),
#        (if (run) mw_ss_code_run(spec) else ""),
        ""), collapse = "\n")
}