diff --git a/R/initOM.R b/R/initOM.R index 6e5111e..5700db3 100644 --- a/R/initOM.R +++ b/R/initOM.R @@ -505,6 +505,116 @@ run_OM <- function(OM_dir, return(dat) } + +#' Initial run of the Base model to get expected and bootstrap values +#' +#' This function is used to initialize the Base model and get expected values +#' for the OM and bootstraps for the EM. +#' @template OM_dir +#' @template verbose +#' @param debug_par_run If set to TRUE, and the run fails, a new folder called +#' error_check will be created, and the model will be run from control start +#' values instead of ss.par. The 2 par files are then compared to help debug +#' the issue with the model run. Defaults to TRUE. +#' @template seed +#' @author Kathryn Doering +#' @importFrom r4ss SS_readdat SS_readstarter SS_writestarter +run_Base <- function(OM_dir, + verbose = FALSE, + debug_par_run = TRUE, + seed = NULL) { + + + if (is.null(seed)) { + seed <- stats::runif(1, 1, 9999999) + } + + start <- r4ss::SS_readstarter(file.path(OM_dir, "starter.ss"), + verbose = FALSE + ) + + start[["init_values_src"]] <- 1 + start[["N_bootstraps"]] <- 2 + start[["seed"]] <- seed + + r4ss::SS_writestarter(start, + dir = OM_dir, verbose = FALSE, overwrite = TRUE + ) + + # run SS and get the data set + run_ss_model(OM_dir, "-maxfn 0 -phase 50 -nohess", + verbose = verbose, + debug_par_run = debug_par_run + ) + + dat <- r4ss::SS_readdat(file.path(OM_dir, "data.ss_new"), + section = 2, + verbose = FALSE + ) + + return(dat) +} + +#' Refit OM model using expected values from original OM +#' +#' This function is used to initialize the Base model and get expected values +#' for the OM and bootstraps for the EM. +#' @template OM_dir +#' @param base_dat Expected Value data list from base model run. +#' @template verbose +#' @param debug_par_run If set to TRUE, and the run fails, a new folder called +#' error_check will be created, and the model will be run from control start +#' values instead of ss.par. The 2 par files are then compared to help debug +#' the issue with the model run. Defaults to TRUE. +#' @template seed +#' @author Kathryn Doering +#' @importFrom r4ss SS_readdat SS_readstarter SS_writestarter +refit_OM <- function(OM_dir, + base_dat, + verbose = FALSE, + debug_par_run = TRUE, + seed = NULL) { + + if (is.null(seed)) { + seed <- stats::runif(1, 1, 9999999) + } + + start <- r4ss::SS_readstarter(file.path(OM_dir, "starter.ss"), + verbose = FALSE + ) + start[["init_values_src"]] <- 1 + start[["detailed_age_structure"]] <- 1 + start[["last_estimation_phase"]] <- 10 + start[["depl_basis"]] <- 1 + start[["depl_denom_frac"]] <- 1 + start[["SPR_basis"]] <- 4 + start[["F_report_units"]] <- 1 + start[["F_report_basis"]] <- 0 + start[["F_age_range"]] <- NULL + start[["ALK_tolerance"]] <- 0 + start[["minyr_sdreport"]] <- -1 + start[["maxyr_sdreport"]] <- -2 + start[["N_bootstraps"]] <- 3 + start[["seed"]] <- seed + + r4ss::SS_writestarter(start, + dir = OM_dir, verbose = FALSE, overwrite = TRUE + ) + + r4ss::SS_writedat(base_dat, file.path(OM_dir, start[["datfile"]]), + overwrite = TRUE, + verbose = FALSE + ) + + # run SS and get the data set + run_ss_model(OM_dir, "-nohess", + verbose = verbose, + debug_par_run = debug_par_run + ) + + return(NULL) +} + #' Get the sampling scheme in a data file. #' #' Determine what the default sampling scheme is for a given data file. diff --git a/R/runSSMSE.R b/R/runSSMSE.R index 5e0bde7..bb34dec 100644 --- a/R/runSSMSE.R +++ b/R/runSSMSE.R @@ -556,6 +556,8 @@ run_SSMSE_iter <- function(out_dir = NULL, OM_in_dir <- out_loc[["OM_in_dir"]] EM_in_dir <- out_loc[["EM_in_dir"]] EM_out_dir <- out_loc[["EM_out_dir"]] + Base_out_dir <- out_loc[["Base_out_dir"]] + Base_in_dir <- out_loc[["Base_in_dir"]] if (!is.null(EM_out_dir)) { EM_out_dir_basename <- strsplit(EM_out_dir, "_init$")[[1]][1] } else { @@ -563,16 +565,32 @@ run_SSMSE_iter <- function(out_dir = NULL, } copy_model_files( OM_in_dir = OM_in_dir, OM_out_dir = OM_out_dir, - EM_in_dir = EM_in_dir, EM_out_dir = EM_out_dir + EM_in_dir = EM_in_dir, EM_out_dir = EM_out_dir, + Base_in_dir = Base_in_dir, Base_out_dir= Base_out_dir ) + # clean model files ---- # want to do this as soon as possible to "fail fast" # for now, this just gets rid of -year value observations, but could do # other things. clean_init_mod_files( OM_out_dir = OM_out_dir, EM_out_dir = EM_out_dir, + Base_out_dir= Base_out_dir, overwrite = TRUE ) + # Add code to create a base model that will be used to sample expected values + # to inform an OM and bootstrapped values to inform the EM's + Base_dat <- run_Base( + OM_dir = Base_out_dir, verbose = verbose, + seed = (iter_seed[["iter"]][1] + 12345) + ) + + browser() + + refit_OM( + OM_dir = OM_out_dir, base_dat = Base_dat, + verbose = verbose, seed = (iter_seed[["iter"]][1] + 12345) + ) # convert sample_struct names ---- # get the full sampling structure for components that the user didnt specify. @@ -588,6 +606,22 @@ run_SSMSE_iter <- function(out_dir = NULL, if(!is.null(sample_struct_hist)){ sample_struct_hist <- convert_to_r4ss_names(sample_struct_hist) } + # Convert the user input parameter modifications into vectors of annual additive deviations + future_base_dat <- convert_future_om_list_to_devs_df(future_om_list = future_om_list, scen_name = scen_name, niter = niter, om_mod_path = Base_out_dir, nyrs = nyrs, global_seed = (iter_seed[["iter"]][1] + 1234)) + + # MSE first iteration ---- + # turn the stock assessment model into an OM + init_base_mod <- create_OM( + OM_out_dir = Base_out_dir, overwrite = TRUE, + sample_struct_hist = sample_struct_hist, + sample_struct = sample_struct, + verbose = verbose, writedat = TRUE, nyrs = nyrs, + nyrs_assess = nyrs_assess, nscen = nscen, + scen_name = scen_name, niter = niter, + future_om_dat = future_base_dat, + seed = (iter_seed[["iter"]][1] + 1234) + ) + # Convert the user input parameter modifications into vectors of annual additive deviations future_om_dat <- convert_future_om_list_to_devs_df(future_om_list = future_om_list, scen_name = scen_name, niter = niter, om_mod_path = OM_out_dir, nyrs = nyrs, global_seed = (iter_seed[["iter"]][1] + 1234)) @@ -606,8 +640,13 @@ run_SSMSE_iter <- function(out_dir = NULL, impl_error <- init_mod[["impl_error"]] # Complete the OM run so it can be use for expect values or bootstrap if (use_SS_boot == TRUE) { + # OM_dat <- run_OM( + # OM_dir = OM_out_dir, boot = use_SS_boot, nboot = 1, + # sample_catch = sample_catch, + # verbose = verbose, init_run = TRUE, seed = (iter_seed[["iter"]][1] + 12345) + # ) OM_dat <- run_OM( - OM_dir = OM_out_dir, boot = use_SS_boot, nboot = 1, + OM_dir = Base_out_dir, boot = use_SS_boot, nboot = 1, sample_catch = sample_catch, verbose = verbose, init_run = TRUE, seed = (iter_seed[["iter"]][1] + 12345) ) diff --git a/R/utils.R b/R/utils.R index 3f55a3c..576511f 100644 --- a/R/utils.R +++ b/R/utils.R @@ -140,11 +140,12 @@ create_scen_list <- function(scen_name_vec, #' #' @template OM_out_dir #' @template EM_out_dir +#' @param Base_out_dir The directory where the base model files will be written. #' @template MS #' @template overwrite #' @importFrom r4ss SS_readstarter SS_readdat SS_writedat -clean_init_mod_files <- function(OM_out_dir, EM_out_dir = NULL, MS = "EM", - overwrite = FALSE) { +clean_init_mod_files <- function(OM_out_dir, EM_out_dir = NULL, Base_out_dir = NULL, + MS = "EM", overwrite = FALSE) { # read in files OM_start <- SS_readstarter(file.path(OM_out_dir, "starter.ss"), verbose = FALSE @@ -157,6 +158,17 @@ clean_init_mod_files <- function(OM_out_dir, EM_out_dir = NULL, MS = "EM", datlist = OM_dat, verbose = FALSE ) + Base_start <- SS_readstarter(file.path(Base_out_dir, "starter.ss"), + verbose = FALSE + ) + Base_dat <- SS_readdat(file.path(Base_out_dir, Base_start[["datfile"]]), + verbose = FALSE, + section = 1 + ) + Base_ctl <- SS_readctl(file.path(Base_out_dir, Base_start[["ctlfile"]]), + datlist = Base_dat, + verbose = FALSE + ) if (!is.null(EM_out_dir)) { EM_start <- SS_readstarter(file.path(EM_out_dir, "starter.ss"), verbose = FALSE @@ -178,9 +190,6 @@ clean_init_mod_files <- function(OM_out_dir, EM_out_dir = NULL, MS = "EM", styr <- OM_dat[["styr"]] endyr <- OM_dat[["endyr"]] - - - # get years in range function get_yrs_in_range <- function(list_name, dat, styr, endyr) { df <- dat[[list_name]] @@ -513,11 +522,17 @@ create_out_dirs <- function(out_dir, niter, OM_name, OM_in_dir, # figure out the OM_name and create the directory within the one just created if (!is.null(OM_name)) { OM_folder_name <- paste0(OM_name, "_OM") + Base_folder_name <- paste0(OM_name, "_OM_Base") } else { OM_folder_name <- paste0(basename(OM_in_dir), "_OM") + Base_folder_name <- paste0(basename(OM_in_dir), "_OM_Base") } OM_out_dir <- file.path(out_dir, OM_folder_name) dir.create(OM_out_dir, showWarnings = FALSE) + + Base_out_dir <- file.path(out_dir, Base_folder_name) + dir.create(Base_out_dir, showWarnings = FALSE) + # Add the EM dir, if necessary if (is.null(EM_name) & !is.null(EM_in_dir)) EM_name <- basename(EM_in_dir) if (!is.null(EM_name) & is.null(EM_in_dir)) { @@ -542,7 +557,8 @@ create_out_dirs <- function(out_dir, niter, OM_name, OM_in_dir, } OM_mod_loc <- list( OM_in_dir = OM_in_dir, OM_out_dir = OM_out_dir, - EM_in_dir = EM_in_dir, EM_out_dir = EM_out_dir + EM_in_dir = EM_in_dir, EM_out_dir = EM_out_dir, + Base_in_dir = OM_in_dir, Base_out_dir = Base_out_dir ) } @@ -587,6 +603,8 @@ locate_in_dirs <- function(OM_name = NULL, OM_in_dir = NULL) { #' @template OM_EM_in_dir #' @template OM_out_dir #' @template EM_out_dir +#' @param Base_in_dir Relative or absolute path to the input base model +#' @param Base_out_dir Relative or absolute path to the base model outputs #' @template verbose #' @return TRUE, if copying is successful #' @@ -594,6 +612,8 @@ copy_model_files <- function(OM_in_dir = NULL, OM_out_dir = NULL, EM_in_dir = NULL, EM_out_dir = NULL, + Base_in_dir = NULL, + Base_out_dir= NULL, verbose = FALSE) { # checks if (!is.null(OM_in_dir)) { @@ -656,7 +676,44 @@ copy_model_files <- function(OM_in_dir = NULL, } else { success_EM <- TRUE } - invisible(c(success_OM = success_OM, success_EM = success_EM)) + + if (!is.null(Base_in_dir)) { + if (!all(c( + "control.ss_new", "data.ss_new", "starter.ss_new", + "forecast.ss_new", "ss.par" + ) %in% list.files(Base_in_dir))) { + stop( + ".ss_new files not found in the original Base directory ", + Base_in_dir, ". Please run the model to make the .ss_new files available." + ) + } + } + # copy over Base ---- + if (!is.null(Base_in_dir) & !is.null(Base_out_dir)) { + if (verbose == TRUE) { + message( + "Copying over .ss_new model files in ", Base_in_dir, + " to ", Base_out_dir, "." + ) + } + success_Base <- r4ss::copy_SS_inputs( + dir.old = Base_in_dir, + dir.new = Base_out_dir, + overwrite = FALSE, + use_ss_new = TRUE, # will rename the ss new files, also. + copy_par = TRUE, + verbose = FALSE + ) + if (success_Base == FALSE) { + stop( + "Problem copying SS Base .ss_new files from ", Base_in_dir, " to ", + Base_out_dir, "." + ) + } + } else { + success_Base <- TRUE + } + invisible(c(success_OM = success_OM, success_EM = success_EM, success_Base = success_Base)) } #' function that creates a combined column to the list_item of interest