From 33a9077e4459cc6b352a1ba5ba7d0f10452b8759 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Mon, 29 Aug 2022 12:04:15 -0700 Subject: [PATCH] refactor: styler and remove OS check --- R/get_settings.R | 51 +++--- R/get_settings_profile.R | 102 ++++++------ R/get_summary.R | 143 +++++++++-------- R/jitter_wrapper.R | 280 ++++++++++++++++++--------------- R/pdf_fxn.R | 19 ++- R/profile_plot.R | 192 ++++++++++++++--------- R/profile_wrapper.R | 325 +++++++++++++++++++------------------- R/rerun_profile_vals.R | 330 ++++++++++++++++++++------------------- R/run_diagnostics.R | 24 +-- 9 files changed, 776 insertions(+), 690 deletions(-) diff --git a/R/get_settings.R b/R/get_settings.R index 2d3482a..15eb553 100644 --- a/R/get_settings.R +++ b/R/get_settings.R @@ -13,15 +13,14 @@ #' #' @examples #' get_settings(list("Njitter" = 10)) -#' +#' #' settings <- list() - +#' get_settings <- function(settings = NULL, verbose = FALSE) { - if (is.vector(settings)) settings <- as.list(settings) - Settings_all = list( - base_name = "model", + Settings_all <- list( + base_name = "model", para_offset = FALSE, run = c("jitter", "profile", "retro"), profile_details = NULL, @@ -29,64 +28,64 @@ get_settings <- function(settings = NULL, verbose = FALSE) { exe = "ss", verbose = FALSE, - # Jitter Settings + # Jitter Settings extras = "-nohess", Njitter = 100, show_in_console = TRUE, printlikes = FALSE, jitter_fraction = 0.05, jitter_init_values_src = NULL, - + # Retrospective Settings oldsubdir = "", - newsubdir = 'retro', + newsubdir = "retro", retro_yrs = -1:-5, overwrite = TRUE, show_in_console = FALSE, - + # Profile Settings remove_files = TRUE, - newctlfile = "control_modified.ss", - linenum = NULL, - string = NULL, + newctlfile = "control_modified.ss", + linenum = NULL, + string = NULL, profilevec = NULL, - usepar = FALSE, - globalpar = FALSE, - parlinenum = NULL, + usepar = FALSE, + globalpar = FALSE, + parlinenum = NULL, parstring = NULL, saveoutput = TRUE, overwrite = TRUE, - whichruns = NULL, + whichruns = NULL, prior_check = FALSE, read_like = TRUE, - init_values_src = 0 + init_values_src = 0 ) - Settings_all$profile_details = get_settings_profile() + Settings_all$profile_details <- get_settings_profile() need <- !names(Settings_all) %in% names(settings) if (verbose) { message("Adding the following objects to settings:\n", paste(names(Settings_all[need]), collapse = "\n"), "\n", - appendLF = TRUE) + appendLF = TRUE + ) } Settings_all <- c(settings, Settings_all[need]) # Check some items - if (!is.null( Settings_all$profile_details) ) { - if (length(Settings_all$profile_details[is.na(Settings_all$profile_details)]) > 0){ + if (!is.null(Settings_all$profile_details)) { + if (length(Settings_all$profile_details[is.na(Settings_all$profile_details)]) > 0) { stop("Missing entry in the get_settings_profile data frame.") } if (!is.numeric(Settings_all$profile_details$low) & - !is.numeric(Settings_all$profile_details$high) & - !is.numeric(Settings_all$profile_details$step_size)){ + !is.numeric(Settings_all$profile_details$high) & + !is.numeric(Settings_all$profile_details$step_size)) { stop("There is a non-numeric value in the low, high, or step size fiedl of the get_settings_profile data frame.") } - if (sum(!Settings_all$profile_details$param_space %in% c("real", "relative", "multiplier")) > 0){ + if (sum(!Settings_all$profile_details$param_space %in% c("real", "relative", "multiplier")) > 0) { stop("The param_space column should be either real or relative in the get_settings_profile data frame.") } - } return(Settings_all) -} \ No newline at end of file +} diff --git a/R/get_settings_profile.R b/R/get_settings_profile.R index 822ea93..8d906fe 100644 --- a/R/get_settings_profile.R +++ b/R/get_settings_profile.R @@ -1,19 +1,19 @@ #' Get Default Settings For Profiles -#' +#' #' Create a matrix of default values for profiling over #' the typical parameters given results will be presented to the #' Pacific Fisheries Management Council. -#' +#' #' The column titled 'param_space' indicated where the range of of the profile parameter #' should be interpretted as relative to the base model estimate vs. across a pre-specified range. #' An example is for R0 where the default setting below indicates that the param_space is relative #' where the low bound for the profile is set = base model log(R0) - 2 and high = base model log(R0) + 2. -#' The default range for M is set as a multiplier to explore a range of (M - 0.40 * M) - (M + 0.40 * M) +#' The default range for M is set as a multiplier to explore a range of (M - 0.40 * M) - (M + 0.40 * M) #' at a step size of 0.005. This range may be too large (or small) with a step size too large (or too small) #' and should be considered if the default settings are appropriate for your specific model. The default setting -#' for steepness is in 'real' space which means that the low and high is in the same parameter space as the -#' parameter. A user can select any of the options for specifying a parameter range for any parameter. -#' +#' for steepness is in 'real' space which means that the low and high is in the same parameter space as the +#' parameter. A user can select any of the options for specifying a parameter range for any parameter. +#' #' @param parameters vector of SS parameter names to conduct a profile for #' @param low a vector of low paramater bounds for the profile #' @param high a vector of upper parameter bounds for the profile @@ -22,9 +22,9 @@ #' real indicates bounds in the parameter space, relative indicates how far to go from the base parameter, and #' multiplier indicates that low and high bounds are set at x\% above and below the base parameter. #' @param use_prior_like options: Option to include or exclude the prior likelihoods in the likelihood -#' profiles. A value of 0 corresponds to the Stock Sythesis were this would exclude -#' and a value of 1 would include the prior likelihood contribution. Parameters with priors used for estimation -#' (e.g., natural mortality, steepness) are often profiled across and including or excluding the prior likelihood +#' profiles. A value of 0 corresponds to the Stock Sythesis were this would exclude +#' and a value of 1 would include the prior likelihood contribution. Parameters with priors used for estimation +#' (e.g., natural mortality, steepness) are often profiled across and including or excluding the prior likelihood #' contribution may be wanted in specific instances. The default setting excludes the prior likelihood contributions. #' #' @return A matrix of low, high, and step size values for the default parameters @@ -34,53 +34,53 @@ #' #' @author Chantel Wetzel & Kelli Johnson #' @export -#' +#' #' @examples #' \dontrun{ -#' +#' #' # Define each parameter in real space -#' get_settings_profile( parameters = c("NatM_p_1_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), -# low = c(0.02, 0.25, 8), -# high = c(0.07, 1.0, 11), -# step_size = c(0.005, 0.05, 0.25), -# param_space = c('real', 'real', 'real'), -#' use_prior_like = c(1, 1, 0)) -#' +#' get_settings_profile( +#' parameters = c("NatM_p_1_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), +#' # low = c(0.02, 0.25, 8), +#' # high = c(0.07, 1.0, 11), +#' # step_size = c(0.005, 0.05, 0.25), +#' # param_space = c('real', 'real', 'real'), +#' use_prior_like = c(1, 1, 0) +#' ) +#' #' # Example 2: Run a profile for natural mortality one with the prior likelihood and one without -#' get_settings_profile( parameters = c("NatM_p_1_Fem_GP_1","NatM_p_1_Fem_GP_1"), -#' low = c(0.40, 0.40), -#' high = c(0.40, 0.40), -#' step_size = c(0.005, 0.005), -#' param_space = c('multiplier', 'multiplier'), -#' use_prior_like = c(0, 1)) -#'} +#' get_settings_profile( +#' parameters = c("NatM_p_1_Fem_GP_1", "NatM_p_1_Fem_GP_1"), +#' low = c(0.40, 0.40), +#' high = c(0.40, 0.40), +#' step_size = c(0.005, 0.005), +#' param_space = c("multiplier", "multiplier"), +#' use_prior_like = c(0, 1) +#' ) +#' } #' -get_settings_profile <- function( parameters = c("NatM_p_1_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), - low = c(0.40, 0.25, -2), - high = c(0.40, 1.0, 2), - step_size = c(0.01, 0.05, 0.25), - param_space = c('multiplier', 'real', 'relative'), - use_prior_like = c(0, 0, 0)) -{ - - if (length(parameters) != length(low) | - length(parameters) != length(high) | - length(parameters) != length(step_size) | - length(parameters) != length(param_space) | - length(parameters) != length(use_prior_like) ){ - stop("Error: input vectors do match in length.") - } +get_settings_profile <- function(parameters = c("NatM_p_1_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), + low = c(0.40, 0.25, -2), + high = c(0.40, 1.0, 2), + step_size = c(0.01, 0.05, 0.25), + param_space = c("multiplier", "real", "relative"), + use_prior_like = c(0, 0, 0)) { + if (length(parameters) != length(low) | + length(parameters) != length(high) | + length(parameters) != length(step_size) | + length(parameters) != length(param_space) | + length(parameters) != length(use_prior_like)) { + stop("Error: input vectors do match in length.") + } - out = data.frame( parameters = parameters, - low = low, - high = high, - step_size = step_size, - param_space = param_space, - use_prior_like = use_prior_like) + out <- data.frame( + parameters = parameters, + low = low, + high = high, + step_size = step_size, + param_space = param_space, + use_prior_like = use_prior_like + ) - return(out) + return(out) } - - - - diff --git a/R/get_summary.R b/R/get_summary.R index f5931cc..44ab6e9 100644 --- a/R/get_summary.R +++ b/R/get_summary.R @@ -1,86 +1,97 @@ #' Generate likelihood profiles #' To be called from the run_diagnostics function after creating #' the model settings using the get_settings function. -#' +#' #' #' @template mydir #' @param name Identify if the csv file show jitter, retro, or profile results -#' @param para SS parameter name that the profile was run across +#' @param para SS parameter name that the profile was run across #' @param vec vector of parameter values the profile covers #' @param profilemodels object created by the SSgetoutput funciton #' @param profilesummary object created by the SSsummarize funciton -#' +#' #' @author Chantel Wetzel & Kelli Johnson #' @export -get_summary <- function(mydir, para, vec, name, profilemodels, profilesummary){ +get_summary <- function(mydir, para, vec, name, profilemodels, profilesummary) { + + # Need to identify a way to determine if a model estmates male growth parameters as offsets from females - # Need to identify a way to determine if a model estmates male growth parameters as offsets from females + # get output + outputs <- profilemodels + quants <- lapply(outputs, "[[", "derived_quants") + status <- sapply(sapply(outputs, "[[", "parameters", simplify = FALSE), "[[", "Status") + bounds <- apply(status, 2, function(x) rownames(outputs[[1]]$parameters)[x %in% c("LO", "HI")]) - # get output - outputs <- profilemodels - quants <- lapply(outputs, "[[", "derived_quants") - status <- sapply(sapply(outputs, "[[", "parameters", simplify = FALSE), "[[", "Status") - bounds <- apply(status, 2, function(x) rownames(outputs[[1]]$parameters)[x %in% c("LO", "HI")]) + out <- data.frame( + "run" = gsub("replist", "", names(outputs)), + "profile_parameter" = para, + "parameter_value" = as.numeric(vec), + "likelihood" = sapply(sapply(outputs, "[[", "likelihoods_used", simplify = FALSE), "[", 1, 1), + "gradient" = sapply(outputs, "[[", "maximum_gradient_component"), + "SB0" = sapply(quants, "[[", "SSB_Virgin", "Value"), + "SBfinal" = sapply(quants, "[[", paste0("SSB_", outputs[[1]]$endyr), "Value"), + "Deplfinal" = sapply(quants, "[[", paste0("Bratio_", outputs[[1]]$endyr), "Value"), + # "Fmsy" = sapply(quants, "[[", "annF_MSY", "Value"), + "Nparsonbounds" = apply(status, 2, function(x) sum(x %in% c("LO", "HI"))), + stringsAsFactors = FALSE + ) - out <- data.frame("run" = gsub("replist", "", names(outputs)), - "profile_parameter" = para, - "parameter_value" = as.numeric(vec), - "likelihood" = sapply(sapply(outputs, "[[", "likelihoods_used", simplify = FALSE), "[", 1, 1), - "gradient" = sapply(outputs, "[[", "maximum_gradient_component"), - "SB0" = sapply(quants, "[[", "SSB_Virgin", "Value"), - "SBfinal" = sapply(quants, "[[", paste0("SSB_", outputs[[1]]$endyr), "Value"), - "Deplfinal" = sapply(quants, "[[", paste0("Bratio_", outputs[[1]]$endyr), "Value"), - #"Fmsy" = sapply(quants, "[[", "annF_MSY", "Value"), - "Nparsonbounds" = apply(status, 2, function(x) sum(x %in% c("LO", "HI"))), - stringsAsFactors = FALSE) + # write tables + write.csv(x = table(unlist(bounds)), file = file.path(mydir, paste0(name, "_parsonbounds.csv")), row.names = FALSE) - # write tables - write.csv(x = table(unlist(bounds)), file= file.path(mydir, paste0(name, "_parsonbounds.csv")), row.names=FALSE) + write.csv(x = out, file = file.path(mydir, paste0(name, "_results.csv")), row.names = FALSE) - write.csv(x = out, file = file.path(mydir, paste0(name, "_results.csv")), row.names = FALSE) + x <- profilesummary + n <- x$n + endyr <- x$endyrs[1] + out <- data.frame( + totlikelihood = as.numeric(x$likelihoods[x$likelihoods$Label == "TOTAL", 1:n]), + surveylike = as.numeric(x$likelihoods[x$likelihoods$Label == "Survey", 1:n]), + discardlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Discard", 1:n]), + lengthlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Length_comp", 1:n]), + agelike = as.numeric(x$likelihoods[x$likelihoods$Label == "Age_comp", 1:n]), + recrlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Recruitment", 1:n]), + forerecrlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Forecast_Recruitment", 1:n]), + priorlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Parm_priors", 1:n]), + parmlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Parm_devs", 1:n]), + R0 = as.numeric(x$pars[x$pars$Label == "SR_LN(R0)", 1:n]), + SB0 = as.numeric(x$SpawnBio[x$SpawnBio$Label == "SSB_Virgin", 1:n]), + SBfinal = as.numeric(x$SpawnBio[x$SpawnBio$Label == paste0("SSB_", endyr), 1:n]), + deplfinal = as.numeric(x$Bratio[x$Bratio$Label == paste0("Bratio_", endyr), 1:n]), + yieldspr = as.numeric(x$quants[x$quants$Label == "Dead_Catch_SPR", 1:n]), + steep = as.numeric(x$pars[x$pars$Label == "SR_BH_steep", 1:n]), + mfem = as.numeric(x$pars[x$pars$Label == "NatM_p_1_Fem_GP_1", 1:n]), + lminfem = as.numeric(x$pars[x$pars$Label == "L_at_Amin_Fem_GP_1", 1:n]), + lmaxfem = as.numeric(x$pars[x$pars$Label == "L_at_Amax_Fem_GP_1", 1:n]), + kfem = as.numeric(x$pars[x$pars$Label == "VonBert_K_Fem_GP_1", 1:n]), + cv1fem = as.numeric(x$pars[grep("young_Fem_GP_1", x$pars$Label), 1:n]), + cv2fem = as.numeric(x$pars[grep("old_Fem_GP_1", x$pars$Label), 1:n]), + mmale = as.numeric(x$pars[x$pars$Label == "NatM_p_1_Mal_GP_1", 1:n]), + lminmale = as.numeric(x$pars[x$pars$Label == "L_at_Amin_Mal_GP_1", 1:n]), + lmaxmale = as.numeric(x$pars[x$pars$Label == "L_at_Amax_Mal_GP_1", 1:n]), + kmale = as.numeric(x$pars[x$pars$Label == "VonBert_K_Mal_GP_1", 1:n]), + cv1male = as.numeric(x$pars[grep("young_Mal_GP_1", x$pars$Label), 1:n]), + cv2male = as.numeric(x$pars[grep("old_Mal_GP_1", x$pars$Label), 1:n]), + stringsAsFactors = FALSE + ) - x <- profilesummary - n <- x$n - endyr <- x$endyrs[1] - out <- data.frame( - totlikelihood = as.numeric(x$likelihoods[x$likelihoods$Label == "TOTAL",1:n]), - surveylike = as.numeric(x$likelihoods[x$likelihoods$Label == "Survey",1:n]), - discardlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Discard",1:n]), - lengthlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Length_comp",1:n]), - agelike = as.numeric(x$likelihoods[x$likelihoods$Label == "Age_comp",1:n]), - recrlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Recruitment",1:n]), - forerecrlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Forecast_Recruitment",1:n]), - priorlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Parm_priors",1:n]), - parmlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Parm_devs",1:n]), - R0 = as.numeric(x$pars[x$pars$Label == "SR_LN(R0)", 1:n]), - SB0 = as.numeric(x$SpawnBio[x$SpawnBio$Label == "SSB_Virgin", 1:n]), - SBfinal = as.numeric(x$SpawnBio[x$SpawnBio$Label == paste0("SSB_", endyr), 1:n]), - deplfinal = as.numeric(x$Bratio[x$Bratio$Label == paste0("Bratio_", endyr), 1:n]), - yieldspr = as.numeric(x$quants[x$quants$Label == "Dead_Catch_SPR", 1:n]), - steep = as.numeric(x$pars[x$pars$Label == "SR_BH_steep", 1:n]), - mfem = as.numeric(x$pars[x$pars$Label == "NatM_p_1_Fem_GP_1", 1:n]), - lminfem = as.numeric(x$pars[x$pars$Label == "L_at_Amin_Fem_GP_1", 1:n]), - lmaxfem = as.numeric(x$pars[x$pars$Label == "L_at_Amax_Fem_GP_1", 1:n]), - kfem = as.numeric(x$pars[x$pars$Label == "VonBert_K_Fem_GP_1", 1:n]), - cv1fem = as.numeric(x$pars[grep("young_Fem_GP_1", x$pars$Label), 1:n]), - cv2fem = as.numeric(x$pars[grep("old_Fem_GP_1", x$pars$Label), 1:n]), - mmale = as.numeric(x$pars[x$pars$Label =="NatM_p_1_Mal_GP_1", 1:n]), - lminmale = as.numeric(x$pars[x$pars$Label =="L_at_Amin_Mal_GP_1", 1:n]), - lmaxmale = as.numeric(x$pars[x$pars$Label =="L_at_Amax_Mal_GP_1", 1:n]), - kmale = as.numeric(x$pars[x$pars$Label =="VonBert_K_Mal_GP_1", 1:n]), - cv1male = as.numeric(x$pars[grep("young_Mal_GP_1", x$pars$Label), 1:n]), - cv2male = as.numeric(x$pars[grep("old_Mal_GP_1", x$pars$Label), 1:n]), - stringsAsFactors = FALSE) + new_out <- t(out) + colnames(new_out) <- vec + if (para == "SR_LN(R0)") { + colnames(new_out) <- paste0("R0 ", vec) + } + if (para == "NatM_p_1_Fem_GP_1") { + colnames(new_out) <- paste0("M_f ", vec) + } + if (para == "NatM_p_1_Mal_GP_1") { + colnames(new_out) <- paste0("M_m ", vec) + } + if (para == "SR_BH_steep") { + colnames(new_out) <- paste0("h ", vec) + } - new_out = t(out) - colnames(new_out) = vec - if (para == "SR_LN(R0)") { colnames(new_out) = paste0("R0 ", vec) } - if (para == "NatM_p_1_Fem_GP_1") { colnames(new_out) = paste0("M_f ", vec) } - if (para == "NatM_p_1_Mal_GP_1") { colnames(new_out) = paste0("M_m ", vec) } - if (para == "SR_BH_steep") { colnames(new_out) = paste0("h ", vec) } - - write.csv(x = new_out, file = file.path(mydir, paste0(name, "_quant_table.csv")), row.names = TRUE) + write.csv(x = new_out, file = file.path(mydir, paste0(name, "_quant_table.csv")), row.names = TRUE) - return () -} \ No newline at end of file + return() +} diff --git a/R/jitter_wrapper.R b/R/jitter_wrapper.R index d006900..d2f97a2 100644 --- a/R/jitter_wrapper.R +++ b/R/jitter_wrapper.R @@ -1,157 +1,178 @@ #' Run [r4ss::jitter] based on `model_settings` -#' +#' #' Code to run jitters for a model #' Output will be saved in an Rdata object called "jitter_output" #' Plots and tables generated to visualize results -#' +#' #' @seealso The following functions interact with `jitter_wrapper`: #' * [run_diagnostics]: calls `jitter_wrapper` #' * [r4ss::jitter]: the workhorse of `jitter_wrapper` that does the jitters #' #' @template mydir #' @template model_settings -#' +#' #' @author Chantel Wetzel #' @return #' Nothing is explicitly returned from `jitter_wrapper`. -#' +#' #' #' @export -jitter_wrapper <- function(mydir, model_settings){ - - if(!file.exists(file.path(mydir, model_settings$base_name, "Report.sso"))) { - message("There is no Report.sso file in the base model directory", file.path(mydir, model_settings$base_name)) - stop() - } - - # Create a jitter folder with the same naming structure as the base model - jitter_dir <- file.path(mydir, paste0(model_settings$base_name, "_jitter_", model_settings$jitter_fraction)) - dir.create(jitter_dir, showWarnings = FALSE) - all_files = list.files(file.path(mydir, model_settings$base_name)) - capture.output(file.copy(from = file.path(mydir, model_settings$base_name, all_files), - to = jitter_dir, - overwrite = TRUE), file = "run_diag_warning.txt") - message("Running jitters: temporarily changing working directory to: ", jitter_dir) - - r4ss::jitter(dir = jitter_dir, - exe = model_settings$exe, - Njitter = model_settings$Njitter, - printlikes = model_settings$printlikes, - verbose = model_settings$verbose, - jitter_fraction = model_settings$jitter_fraction, - init_values_src = model_settings$jitter_init_values_src, - extras = model_settings$extras, - show_in_console = model_settings$show_in_console) +jitter_wrapper <- function(mydir, model_settings) { + if (!file.exists(file.path(mydir, model_settings$base_name, "Report.sso"))) { + message("There is no Report.sso file in the base model directory", file.path(mydir, model_settings$base_name)) + stop() + } + + # Create a jitter folder with the same naming structure as the base model + jitter_dir <- file.path(mydir, paste0(model_settings$base_name, "_jitter_", model_settings$jitter_fraction)) + dir.create(jitter_dir, showWarnings = FALSE) + all_files <- list.files(file.path(mydir, model_settings$base_name)) + capture.output(file.copy( + from = file.path(mydir, model_settings$base_name, all_files), + to = jitter_dir, + overwrite = TRUE + ), file = "run_diag_warning.txt") + message("Running jitters: temporarily changing working directory to: ", jitter_dir) + + r4ss::jitter( + dir = jitter_dir, + exe = model_settings$exe, + Njitter = model_settings$Njitter, + printlikes = model_settings$printlikes, + verbose = model_settings$verbose, + jitter_fraction = model_settings$jitter_fraction, + init_values_src = model_settings$jitter_init_values_src, + extras = model_settings$extras, + show_in_console = model_settings$show_in_console + ) #### Read in results using other r4ss functions - #keys <- gsub("Report([0-9]+)\\.sso", "\\1", dir(jitter_dir, pattern = "^Report[0-9]")) - #keys <- type.convert(keys)[order(type.convert(keys))] - keys <- 1:model_settings$Njitter - profilemodels <- r4ss::SSgetoutput(dirvec = jitter_dir, - keyvec = keys, - getcovar = FALSE, - forecast = FALSE, - verbose = FALSE, - listlists = TRUE, - underscore = FALSE, - save.lists = FALSE) - - # summarize output - profilesummary <- r4ss::SSsummarize(profilemodels) - - # Read in the base model - base <- r4ss::SS_output(dir = file.path(mydir, model_settings$base_name), - covar = FALSE, - printstats = FALSE, - verbose = FALSE) + # keys <- gsub("Report([0-9]+)\\.sso", "\\1", dir(jitter_dir, pattern = "^Report[0-9]")) + # keys <- type.convert(keys)[order(type.convert(keys))] + keys <- 1:model_settings$Njitter + profilemodels <- r4ss::SSgetoutput( + dirvec = jitter_dir, + keyvec = keys, + getcovar = FALSE, + forecast = FALSE, + verbose = FALSE, + listlists = TRUE, + underscore = FALSE, + save.lists = FALSE + ) - est <- base$likelihoods_used[1, 1] - like <- as.numeric(profilesummary$likelihoods[1, keys]) - ymax <- as.numeric(quantile(unlist(profilesummary$likelihoods[1, keys]), 0.80)) - ymin <- min(like - est) + 1 + # summarize output + profilesummary <- r4ss::SSsummarize(profilemodels) - jitter_output <- list() - jitter_output$plotdir <- jitter_dir - jitter_output$est <- est - jitter_output$keys <- keys - jitter_output$like <- like - jitter_output$model_settings <- model_settings - jitter_output$profilesummary <- profilesummary - jitter_output$profilemodels <- profilemodels + # Read in the base model + base <- r4ss::SS_output( + dir = file.path(mydir, model_settings$base_name), + covar = FALSE, + printstats = FALSE, + verbose = FALSE + ) - save( - jitter_dir, - est, - keys, - like, - model_settings, - profilesummary, - profilemodels, - file = file.path(jitter_dir, "jitter_output.Rdata") - ) + est <- base$likelihoods_used[1, 1] + like <- as.numeric(profilesummary$likelihoods[1, keys]) + ymax <- as.numeric(quantile(unlist(profilesummary$likelihoods[1, keys]), 0.80)) + ymin <- min(like - est) + 1 + + jitter_output <- list() + jitter_output$plotdir <- jitter_dir + jitter_output$est <- est + jitter_output$keys <- keys + jitter_output$like <- like + jitter_output$model_settings <- model_settings + jitter_output$profilesummary <- profilesummary + jitter_output$profilemodels <- profilemodels + + save( + jitter_dir, + est, + keys, + like, + model_settings, + profilesummary, + profilemodels, + file = file.path(jitter_dir, "jitter_output.Rdata") + ) ylab <- "Change in negative log-likelihood" xlab <- "Iteration" - pngfun(wd = jitter_dir, file = paste0("Jitter_", model_settings$jitter_fraction, '.png'), h = 12, w = 9) - plot(keys, like - est, ylim = c(ymin, ymax), cex.axis = 1.25, cex.lab = 1.25, - ylab = ylab, xlab = xlab) - abline(h = 0, col = 'darkgrey', lwd = 2) - find = which(est == like) - points(keys[find], (like - est)[find], col = 'green3', pch = 16, cex = 1.1) - find = which(like - est > 0) - points(keys[find], (like - est)[find], col = 'blue', pch = 16) - if (sum(like - est < 0) > 0) { - find = like - est < 0 - points(keys[find], (like - est)[find], col = 'red', pch = 16, cex = 1.1) - mtext(side = 3, cex = 1.25, - "Warning: A lower NLL was found. Update your base model.") - } - legend('topleft', legend = c("Base Model Likelihood", "Higher Likelihood", "Lower Likelihood"), - bty = 'n', pch = 16, col = c('green3', 'blue', 'red')) - dev.off() - - if (ymax > 100){ - pngfun(wd = jitter_dir, file = paste0("Jitter_Zoomed_SubPlot_", model_settings$jitter_fraction, '.png'), h = 12, w = 9) - plot(keys, like - est, ylim = c(ymin, 100), cex.axis = 1.25, cex.lab = 1.25, - ylab = ylab, xlab = xlab) - abline(h = 0, col = 'darkgrey', lwd = 2) - find = which(est == like) - points(keys[find], (like - est)[find], col = 'green3', pch = 16, cex = 1.1) - find = which(like - est > 0) - points(keys[find], (like - est)[find], col = 'blue', pch = 16) - if (sum(like - est < 0) > 0) { - find = like - est < 0 - points(keys[find], (like - est)[find], col = 'red', pch = 16, cex = 1.1) - mtext(side = 3, cex = 1.25, - "Warning: Only jitters near the base model shown") - } - legend('topleft', legend = c("Base Model Likelihood", "Higher Likelihood", "Lower Likelihood"), - bty = 'n', pch = 16, col = c('green3', 'blue', 'red')) - dev.off() - - } - - # get output - outputs <- profilemodels - quants <- lapply(outputs, "[[", "derived_quants") - status <- sapply(sapply(outputs, "[[", "parameters", simplify = FALSE), "[[", "Status") - bounds <- apply(status, 2, function(x) rownames(outputs[[1]]$parameters)[x %in% c("LO", "HI")]) - out <- data.frame("run" = gsub("replist", "", names(outputs)), - "likelihood" = sapply(sapply(outputs, "[[", "likelihoods_used", simplify = FALSE), "[", 1, 1), - "gradient" = sapply(outputs, "[[", "maximum_gradient_component"), - "SB0" = sapply(quants, "[[", "SSB_Virgin", "Value"), - "SBfinal" = sapply(quants, "[[", paste0("SSB_", profilesummary$endyrs[1]), "Value"), - "Nparsonbounds" = apply(status, 2, function(x) sum(x %in% c("LO", "HI"))), - "Lowest NLL" = ifelse(min(like) == like, "Best Fit", 0), - stringsAsFactors = FALSE) + pngfun(wd = jitter_dir, file = paste0("Jitter_", model_settings$jitter_fraction, ".png"), h = 12, w = 9) + plot(keys, like - est, + ylim = c(ymin, ymax), cex.axis = 1.25, cex.lab = 1.25, + ylab = ylab, xlab = xlab + ) + abline(h = 0, col = "darkgrey", lwd = 2) + find <- which(est == like) + points(keys[find], (like - est)[find], col = "green3", pch = 16, cex = 1.1) + find <- which(like - est > 0) + points(keys[find], (like - est)[find], col = "blue", pch = 16) + if (sum(like - est < 0) > 0) { + find <- like - est < 0 + points(keys[find], (like - est)[find], col = "red", pch = 16, cex = 1.1) + mtext( + side = 3, cex = 1.25, + "Warning: A lower NLL was found. Update your base model." + ) + } + legend("topleft", + legend = c("Base Model Likelihood", "Higher Likelihood", "Lower Likelihood"), + bty = "n", pch = 16, col = c("green3", "blue", "red") + ) + dev.off() + + if (ymax > 100) { + pngfun(wd = jitter_dir, file = paste0("Jitter_Zoomed_SubPlot_", model_settings$jitter_fraction, ".png"), h = 12, w = 9) + plot(keys, like - est, + ylim = c(ymin, 100), cex.axis = 1.25, cex.lab = 1.25, + ylab = ylab, xlab = xlab + ) + abline(h = 0, col = "darkgrey", lwd = 2) + find <- which(est == like) + points(keys[find], (like - est)[find], col = "green3", pch = 16, cex = 1.1) + find <- which(like - est > 0) + points(keys[find], (like - est)[find], col = "blue", pch = 16) + if (sum(like - est < 0) > 0) { + find <- like - est < 0 + points(keys[find], (like - est)[find], col = "red", pch = 16, cex = 1.1) + mtext( + side = 3, cex = 1.25, + "Warning: Only jitters near the base model shown" + ) + } + legend("topleft", + legend = c("Base Model Likelihood", "Higher Likelihood", "Lower Likelihood"), + bty = "n", pch = 16, col = c("green3", "blue", "red") + ) + dev.off() + } + + # get output + outputs <- profilemodels + quants <- lapply(outputs, "[[", "derived_quants") + status <- sapply(sapply(outputs, "[[", "parameters", simplify = FALSE), "[[", "Status") + bounds <- apply(status, 2, function(x) rownames(outputs[[1]]$parameters)[x %in% c("LO", "HI")]) + out <- data.frame( + "run" = gsub("replist", "", names(outputs)), + "likelihood" = sapply(sapply(outputs, "[[", "likelihoods_used", simplify = FALSE), "[", 1, 1), + "gradient" = sapply(outputs, "[[", "maximum_gradient_component"), + "SB0" = sapply(quants, "[[", "SSB_Virgin", "Value"), + "SBfinal" = sapply(quants, "[[", paste0("SSB_", profilesummary$endyrs[1]), "Value"), + "Nparsonbounds" = apply(status, 2, function(x) sum(x %in% c("LO", "HI"))), + "Lowest NLL" = ifelse(min(like) == like, "Best Fit", 0), + stringsAsFactors = FALSE + ) # Write a md file to be included in a stock assessment document # Text was pirated from @chantelwetzel-noaa's 2021 dover assessment file_md <- file.path(jitter_dir, "model-results-jitter.md") sink(file_md) on.exit(sink(), add = TRUE) - cat(sep = "", + cat( + sep = "", "Model convergence was in part based on starting the minimization process ", "from dispersed values of the maximum likelihood estimates to determine if the ", "estimation routine results in a smaller likelihood.\n", @@ -175,10 +196,9 @@ jitter_wrapper <- function(mydir, model_settings){ "best fit to the data given the assumptions made.\n" ) - # write tables - write.csv(x = table(unlist(bounds)), file = file.path(jitter_dir, "jitter_parsonbounds.csv"), row.names = FALSE) - write.csv(x = out, file = file.path(jitter_dir, "jitter_results.csv"), row.names = FALSE) + # write tables + write.csv(x = table(unlist(bounds)), file = file.path(jitter_dir, "jitter_parsonbounds.csv"), row.names = FALSE) + write.csv(x = out, file = file.path(jitter_dir, "jitter_results.csv"), row.names = FALSE) - message("Finished jitters.") - -} \ No newline at end of file + message("Finished jitters.") +} diff --git a/R/pdf_fxn.R b/R/pdf_fxn.R index 905479d..75bda64 100644 --- a/R/pdf_fxn.R +++ b/R/pdf_fxn.R @@ -8,15 +8,14 @@ #' @author by Chantel Wetzel #' @export pngfun <- function(wd, file, w = 7, h = 7, pt = 12) { - file <- file.path(wd, file) - cat('writing PNG to',file,'\n') - png(filename = file, - width = w, - height = h, - units = 'in', - res = 300, - pointsize = pt + cat("writing PNG to", file, "\n") + png( + filename = file, + width = w, + height = h, + units = "in", + res = 300, + pointsize = pt ) - -} \ No newline at end of file +} diff --git a/R/profile_plot.R b/R/profile_plot.R index 2513ed1..2901647 100644 --- a/R/profile_plot.R +++ b/R/profile_plot.R @@ -17,81 +17,115 @@ #' @export #' @seealso [profile_wrapper] and [rerun_profile_vals] call `profile_plot`. -profile_plot <- function(mydir, rep, para, profilesummary){ - +profile_plot <- function(mydir, rep, para, profilesummary) { label <- ifelse(para == "SR_LN(R0)", expression(log(italic(R)[0])), - ifelse(para %in% c("NatM_p_1_Fem_GP_1", "NatM_uniform_Fem_GP_1"), "Natural Mortality (female)", - ifelse(para %in% c("NatM_p_1_Mal_GP_1", "NatM_uniform_Mal_GP_1"), "Natural Mortality (male)", - ifelse(para == "SR_BH_steep", expression(Steepness~(italic(h))), - para)))) + ifelse(para %in% c("NatM_p_1_Fem_GP_1", "NatM_uniform_Fem_GP_1"), "Natural Mortality (female)", + ifelse(para %in% c("NatM_p_1_Mal_GP_1", "NatM_uniform_Mal_GP_1"), "Natural Mortality (male)", + ifelse(para == "SR_BH_steep", expression(Steepness ~ (italic(h))), + para + ) + ) + ) + ) - get = ifelse(para == "SR_LN(R0)", "R0", para) + get <- ifelse(para == "SR_LN(R0)", "R0", para) - if(para %in% c("SR_LN(R0)", "NatM_p_1_Fem_GP_1", "NatM_p_1_Mal_GP_1", - "NatM_uniform_Fem_GP_1", "NatM_uniform_Mal_GP_1", "SR_BH_steep")) { - exact = FALSE + if (para %in% c( + "SR_LN(R0)", "NatM_p_1_Fem_GP_1", "NatM_p_1_Mal_GP_1", + "NatM_uniform_Fem_GP_1", "NatM_uniform_Mal_GP_1", "SR_BH_steep" + )) { + exact <- FALSE } else { - exact = TRUE - } + exact <- TRUE + } n <- 1:profilesummary$n - - like_comp <- unique(profilesummary$likelihoods_by_fleet$Label[ - c(-grep("_lambda", profilesummary$likelihoods_by_fleet$Label), - -grep("_N_use", profilesummary$likelihoods_by_fleet$Label), - -grep("_N_skip", profilesummary$likelihoods_by_fleet$Label))]) + + like_comp <- unique(profilesummary$likelihoods_by_fleet$Label[ + c( + -grep("_lambda", profilesummary$likelihoods_by_fleet$Label), + -grep("_N_use", profilesummary$likelihoods_by_fleet$Label), + -grep("_N_skip", profilesummary$likelihoods_by_fleet$Label) + ) + ]) ii <- which(profilesummary$likelihoods_by_fleet$Label %in% like_comp) - check <- aggregate(ALL~Label, profilesummary$likelihoods_by_fleet[ii,], FUN = sum) - use <- check[which(check$ALL != 0),"Label"] + check <- aggregate(ALL ~ Label, profilesummary$likelihoods_by_fleet[ii, ], FUN = sum) + use <- check[which(check$ALL != 0), "Label"] # If present remove the likes that we don't typically show use <- use[which(!use %in% c("Disc_like", "Catch_like", "mnwt_like"))] - tot_plot <- length(use) - if(tot_plot == 1) { panel <- c(2, 1) } - if(tot_plot != 1 & tot_plot <= 3) { panel <- c(3, 1) } - if(tot_plot >= 3) { panel <- c(2, 2) } - if(tot_plot >= 4) { panel <- c(3, 2) } + tot_plot <- length(use) + if (tot_plot == 1) { + panel <- c(2, 1) + } + if (tot_plot != 1 & tot_plot <= 3) { + panel <- c(3, 1) + } + if (tot_plot >= 3) { + panel <- c(2, 2) + } + if (tot_plot >= 4) { + panel <- c(3, 2) + } # Determine the y-axis for the profile plot for all data types together ymax1 <- max(profilesummary$likelihoods[1, n]) - min(profilesummary$likelihoods[1, n]) - if(ymax1 > 70) { ymax1 = 70} - if(ymax1 < 5) { ymax1 = 5} + if (ymax1 > 70) { + ymax1 <- 70 + } + if (ymax1 < 5) { + ymax1 <- 5 + } # Determine the y-axis for the piner profile plots by each data type lab.row <- ncol(profilesummary$likelihoods) - ymax2 <- max(apply(profilesummary$likelihoods[-1, -lab.row], 1, max) - - apply(profilesummary$likelihoods[-1, -lab.row], 1, min)) - if(ymax2 > 70) { ymax2 = 70} - if(ymax2 < 5) { ymax2 = 5} + ymax2 <- max(apply(profilesummary$likelihoods[-1, -lab.row], 1, max) - + apply(profilesummary$likelihoods[-1, -lab.row], 1, min)) + if (ymax2 > 70) { + ymax2 <- 70 + } + if (ymax2 < 5) { + ymax2 <- 5 + } - pngfun(wd = mydir, file = paste0("piner_panel_", para, ".png"), h= 7, w = 7) + pngfun(wd = mydir, file = paste0("piner_panel_", para, ".png"), h = 7, w = 7) par(mfrow = panel) - r4ss::SSplotProfile(summaryoutput = profilesummary, main = "Changes in total likelihood", profile.string = get, - profile.label = label, ymax = ymax1, exact = exact) - abline(h = 1.92, lty = 3, col = 'red') - - if ("Length_like" %in% use){ - r4ss::PinerPlot (summaryoutput = profilesummary, plot = TRUE, print = FALSE, component = "Length_like", - main = "Length-composition likelihoods", profile.string = get, profile.label = label, - exact = exact, ylab = "Change in -log-likelihood", legendloc = "topright", ymax = ymax2) + r4ss::SSplotProfile( + summaryoutput = profilesummary, main = "Changes in total likelihood", profile.string = get, + profile.label = label, ymax = ymax1, exact = exact + ) + abline(h = 1.92, lty = 3, col = "red") + + if ("Length_like" %in% use) { + r4ss::PinerPlot( + summaryoutput = profilesummary, plot = TRUE, print = FALSE, component = "Length_like", + main = "Length-composition likelihoods", profile.string = get, profile.label = label, + exact = exact, ylab = "Change in -log-likelihood", legendloc = "topright", ymax = ymax2 + ) } - - if ("Age_like" %in% use){ - r4ss::PinerPlot (summaryoutput = profilesummary, plot = TRUE, print = FALSE, component = "Age_like", - main = "Age-composition likelihoods", profile.string = get, profile.label = label, - exact = exact, ylab = "Change in -log-likelihood", legendloc = "topright", ymax = ymax2) + + if ("Age_like" %in% use) { + r4ss::PinerPlot( + summaryoutput = profilesummary, plot = TRUE, print = FALSE, component = "Age_like", + main = "Age-composition likelihoods", profile.string = get, profile.label = label, + exact = exact, ylab = "Change in -log-likelihood", legendloc = "topright", ymax = ymax2 + ) } - - if ("Surv_like" %in% use){ - r4ss::PinerPlot (summaryoutput = profilesummary, plot = TRUE, print = FALSE, component = "Surv_like", - main = "Survey likelihoods", profile.string = get, profile.label = label, - exact = exact, ylab = "Change in -log-likelihood", legendloc = "topright", ymax = ymax2) + + if ("Surv_like" %in% use) { + r4ss::PinerPlot( + summaryoutput = profilesummary, plot = TRUE, print = FALSE, component = "Surv_like", + main = "Survey likelihoods", profile.string = get, profile.label = label, + exact = exact, ylab = "Change in -log-likelihood", legendloc = "topright", ymax = ymax2 + ) } - if ("Init_equ_like" %in% use){ - r4ss::PinerPlot (summaryoutput = profilesummary, plot = TRUE, print = FALSE, component = "Init_equ_like", - main = "Initial equilibrium likelihoods", profile.string = get, profile.label = label, - exact = exact, ylab = "Change in -log-likelihood", legendloc = "topright", ymax = ymax2) + if ("Init_equ_like" %in% use) { + r4ss::PinerPlot( + summaryoutput = profilesummary, plot = TRUE, print = FALSE, component = "Init_equ_like", + main = "Initial equilibrium likelihoods", profile.string = get, profile.label = label, + exact = exact, ylab = "Change in -log-likelihood", legendloc = "topright", ymax = ymax2 + ) } dev.off() @@ -99,38 +133,41 @@ profile_plot <- function(mydir, rep, para, profilesummary){ maxyr <- min(profilesummary$endyrs) minyr <- max(profilesummary$startyrs) - est <- rep$parameters[rep$parameters$Label == para, "Value", 2] - sb0_est <- rep$derived_quants[rep$derived_quants$Label == "SSB_Virgin", "Value"] - sbf_est <- rep$derived_quants[rep$derived_quants$Label == paste0("SSB_", maxyr), "Value"] + est <- rep$parameters[rep$parameters$Label == para, "Value", 2] + sb0_est <- rep$derived_quants[rep$derived_quants$Label == "SSB_Virgin", "Value"] + sbf_est <- rep$derived_quants[rep$derived_quants$Label == paste0("SSB_", maxyr), "Value"] depl_est <- rep$derived_quants[rep$derived_quants$Label == paste0("Bratio_", maxyr), "Value"] - x <- as.numeric(profilesummary$pars[profilesummary$pars$Label == para, n]) - like <- as.numeric(profilesummary$likelihoods[profilesummary$likelihoods$Label == "TOTAL", n] - rep$likelihoods_used[1,1]) - ylike<- c(min(like) + ifelse(min(like) != 0, -0.5, 0), max(like)) - sb0 <- as.numeric(profilesummary$SpawnBio[na.omit(profilesummary$SpawnBio$Label) == "SSB_Virgin", n]) - sbf <- as.numeric(profilesummary$SpawnBio[na.omit(profilesummary$SpawnBio$Yr) == maxyr, n]) + x <- as.numeric(profilesummary$pars[profilesummary$pars$Label == para, n]) + like <- as.numeric(profilesummary$likelihoods[profilesummary$likelihoods$Label == "TOTAL", n] - rep$likelihoods_used[1, 1]) + ylike <- c(min(like) + ifelse(min(like) != 0, -0.5, 0), max(like)) + sb0 <- as.numeric(profilesummary$SpawnBio[na.omit(profilesummary$SpawnBio$Label) == "SSB_Virgin", n]) + sbf <- as.numeric(profilesummary$SpawnBio[na.omit(profilesummary$SpawnBio$Yr) == maxyr, n]) depl <- as.numeric(profilesummary$Bratio[na.omit(profilesummary$Bratio$Yr) == maxyr, n]) - # Get the relative management targets - only grab the first element since the targets should be the same - btarg <- as.numeric(profilesummary$btargs[1]) + # Get the relative management targets - only grab the first element since the targets should be the same + btarg <- as.numeric(profilesummary$btargs[1]) thresh <- ifelse(btarg == 0.40, 0.25, 0.125) pngfun(wd = mydir, file = paste0("parameter_panel_", para, ".png"), h = 7, w = 7) - par(mfrow = c(2,2), mar = c(4,4,2,2), oma = c(1,1,1,1)) + par(mfrow = c(2, 2), mar = c(4, 4, 2, 2), oma = c(1, 1, 1, 1)) # parameter vs. likelihood plot(x, like, type = "l", lwd = 2, xlab = label, ylab = "Change in -log-likelihood", ylim = ylike) - abline(h = 0, lty = 2, col = 'black') - if(max(ylike) < 40) { - abline(h = 1.92, lty = 3, col = 'red') - abline(h = -1.92, lty = 3, col = 'red') } + abline(h = 0, lty = 2, col = "black") + if (max(ylike) < 40) { + abline(h = 1.92, lty = 3, col = "red") + abline(h = -1.92, lty = 3, col = "red") + } points(est, 0, pch = 21, col = "black", bg = "blue", cex = 1.5) # parameter vs. final depletion plot(x, depl, type = "l", lwd = 2, xlab = label, ylab = "Fraction of unfished", ylim = c(0, 1.2)) points(est, depl_est, pch = 21, col = "black", bg = "blue", cex = 1.5) abline(h = c(thresh, btarg), lty = c(2, 2), col = c("red", "darkgreen")) - legend("bottomright", legend = c("Management target", "Minimum stock size threshold"), - lty = 2, col = c('red', 'darkgreen'), bty = 'n') + legend("bottomright", + legend = c("Management target", "Minimum stock size threshold"), + lty = 2, col = c("red", "darkgreen"), bty = "n" + ) # parameter vs. SB0 plot(x, sb0, type = "l", lwd = 2, xlab = label, ylab = expression(SB[0]), ylim = c(0, max(sb0))) @@ -144,11 +181,15 @@ profile_plot <- function(mydir, rep, para, profilesummary){ # Create the sb and depl trajectories plot # Figure out what the base model parameter is in order to label that in the plot - get = ifelse(para == "SR_LN(R0)", "log(R0)", - ifelse(para %in% c("NatM_uniform_Fem_GP_1", "NatM_p_1_Fem_GP_1"), "M (f)", - ifelse(para %in% c("NatM_uniform_Mal_GP_1", "NatM_p_1_Mal_GP_1"), "M (m)", - ifelse(para == "SR_BH_steep", "h", - para)))) + get <- ifelse(para == "SR_LN(R0)", "log(R0)", + ifelse(para %in% c("NatM_uniform_Fem_GP_1", "NatM_p_1_Fem_GP_1"), "M (f)", + ifelse(para %in% c("NatM_uniform_Mal_GP_1", "NatM_p_1_Mal_GP_1"), "M (m)", + ifelse(para == "SR_BH_steep", "h", + para + ) + ) + ) + ) r4ss::SSplotComparisons( summaryoutput = profilesummary, @@ -163,5 +204,4 @@ profile_plot <- function(mydir, rep, para, profilesummary){ pdf = FALSE, print = TRUE, plot = FALSE, filenameprefix = paste0(para, "_trajectories_") ) - } diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index 1a43d6b..a795826 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -3,7 +3,7 @@ #' Generate likelihood profiles #' To be called from the run_diagnostics function after creating #' the model settings using the get_settings function. -#' +#' #' #' @seealso The following functions interact with `profile_wrapper`: #' * [run_diagnostics]: calls `profile_wrapper` @@ -12,201 +12,210 @@ #' #' @template mydir #' @template model_settings -#' +#' #' @author Chantel Wetzel. #' @return #' Nothing is explicitly returned from `profile_wrapper`. #' The following objects are saved to the disk. #' @export -profile_wrapper <- function(mydir, model_settings){ +profile_wrapper <- function(mydir, model_settings) { # Add the round_any function from the plyr package to avoid conflicts between # plyr and dplyr. - round_any <- function(x, accuracy, f = round){ - f(x / accuracy) * accuracy} - - # figure out name of executable based on 'exe' input which may contain .exe - if(length(grep(".exe", tolower(file.path(mydir, model_settings$base_name)))) == 1) { - # if input 'exe' includes .exe then assume it's Windows and just use the name - exe <- model_settings$exe - } else { - # if 'model' doesn't include .exe then append it (for Windows computers only) - exe <- paste(model_settings$exe, ifelse(OS == "Windows", ".exe", ""), sep="") + round_any <- function(x, accuracy, f = round) { + f(x / accuracy) * accuracy } + + check_exe <- paste0(model_settings$exe, ".exe") # check whether exe is in directory - if(!exe %in% list.files(file.path(mydir, model_settings$base_name))) { - stop("Executable ", exe, " not found in ", file.path(mydir, model_settings$base_name)) + if (!check_exe %in% list.files(file.path(mydir, model_settings$base_name))) { + stop("Executable not found in ", file.path(mydir, model_settings$base_name)) } N <- nrow(model_settings$profile_details) - for (aa in 1:N){ - + for (aa in 1:N) { para <- model_settings$profile_details$parameters[aa] prior_used <- model_settings$profile_details$use_prior_like[aa] # Create a profile folder with the same naming structure as the base model # Add a label to show if prior was used or not - profile_dir <- file.path(mydir, paste0(model_settings$base_name, "_profile_", para, "_prior_like_", prior_used)) - dir.create(profile_dir, showWarnings = FALSE) + profile_dir <- file.path(mydir, paste0(model_settings$base_name, "_profile_", para, "_prior_like_", prior_used)) + dir.create(profile_dir, showWarnings = FALSE) - # Check for existing files and delete - if (model_settings$remove_files & length(list.files(profile_dir)) != 0) { - remove <- list.files(profile_dir) - file.remove(file.path(profile_dir, remove)) - } + # Check for existing files and delete + if (model_settings$remove_files & length(list.files(profile_dir)) != 0) { + remove <- list.files(profile_dir) + file.remove(file.path(profile_dir, remove)) + } - all_files <- list.files(file.path(mydir, model_settings$base_name)) - capture.output(file.copy(from = file.path(mydir, model_settings$base_name, all_files), - to = profile_dir, overwrite = TRUE), file = "run_diag_warning.txt") - message(paste0( "Running profile for ", para, ".") ) + all_files <- list.files(file.path(mydir, model_settings$base_name)) + capture.output(file.copy( + from = file.path(mydir, model_settings$base_name, all_files), + to = profile_dir, overwrite = TRUE + ), file = "run_diag_warning.txt") + message(paste0("Running profile for ", para, ".")) - # Copy the control file to run from the copy + # Copy the control file to run from the copy if (!file.exists(file.path(profile_dir, "control.ss_new"))) { orig_dir <- getwd() setwd(profile_dir) - r4ss::run(dir = profile_dir, - exe = model_settings$exe, - extras = model_settings$extras) + r4ss::run( + dir = profile_dir, + exe = model_settings$exe, + extras = model_settings$extras + ) setwd(orig_dir) } - # Use the SS_parlines funtion to ensure that the input parameter can be found - check_para <- r4ss::SS_parlines( - ctlfile = "control.ss_new", - dir = profile_dir, - verbose = FALSE, - version = model_settings$version, - active = FALSE)$Label == para + # Use the SS_parlines funtion to ensure that the input parameter can be found + check_para <- r4ss::SS_parlines( + ctlfile = "control.ss_new", + dir = profile_dir, + verbose = FALSE, + version = model_settings$version, + active = FALSE + )$Label == para - if(sum(check_para) == 0) { + if (sum(check_para) == 0) { print(para) - stop(paste0( "The input profile_custom does not match a parameter in the control.ss_new file.")) - } - + stop(paste0("The input profile_custom does not match a parameter in the control.ss_new file.")) + } + file.copy(file.path(profile_dir, "control.ss_new"), file.path(profile_dir, model_settings$newctlfile)) # Change the control file name in the starter file - starter <- r4ss::SS_readstarter(file = file.path(profile_dir, 'starter.ss')) - starter$ctlfile <- "control_modified.ss" - starter$init_values_src <- model_settings$init_values_src - # make sure the prior likelihood is calculated for non-estimated quantities - starter$prior_like <- prior_used - r4ss::SS_writestarter(mylist = starter, dir = profile_dir, overwrite = TRUE) + starter <- r4ss::SS_readstarter(file = file.path(profile_dir, "starter.ss")) + starter$ctlfile <- "control_modified.ss" + starter$init_values_src <- model_settings$init_values_src + # make sure the prior likelihood is calculated for non-estimated quantities + starter$prior_like <- prior_used + r4ss::SS_writestarter(mylist = starter, dir = profile_dir, overwrite = TRUE) # Read in the base model - rep <- r4ss::SS_output(file.path(mydir, model_settings$base_name), covar = FALSE, - printstats = FALSE, verbose = FALSE) - est <- rep$parameters[rep$parameters$Label == para, "Value"] - - # Determine the parameter range - if (model_settings$profile_details$param_space[aa] == 'relative') { - range <- c( est + model_settings$profile_details$low[aa], - est + model_settings$profile_details$high[aa] ) - } - if (model_settings$profile_details$param_space[aa] == 'multiplier') { - range <- c( est - est * model_settings$profile_details$low[aa], - est + est * model_settings$profile_details$high[aa] ) + rep <- r4ss::SS_output(file.path(mydir, model_settings$base_name), + covar = FALSE, + printstats = FALSE, verbose = FALSE + ) + est <- rep$parameters[rep$parameters$Label == para, "Value"] + + # Determine the parameter range + if (model_settings$profile_details$param_space[aa] == "relative") { + range <- c( + est + model_settings$profile_details$low[aa], + est + model_settings$profile_details$high[aa] + ) + } + if (model_settings$profile_details$param_space[aa] == "multiplier") { + range <- c( + est - est * model_settings$profile_details$low[aa], + est + est * model_settings$profile_details$high[aa] + ) + } + if (model_settings$profile_details$param_space[aa] == "real") { + range <- c( + model_settings$profile_details$low[aa], + model_settings$profile_details$high[aa] + ) } - if (model_settings$profile_details$param_space[aa] == 'real') { - range <- c(model_settings$profile_details$low[aa], - model_settings$profile_details$high[aa]) - } - step_size <- model_settings$profile_details$step_size[aa] - - # Create parameter vect from base down and the base up - if (est != round_any(est, step_size, f = floor)) { - low <- rev(seq(round_any(range[1], step_size, f = ceiling), - round_any(est, step_size, f = floor), step_size)) + step_size <- model_settings$profile_details$step_size[aa] + + # Create parameter vect from base down and the base up + if (est != round_any(est, step_size, f = floor)) { + low <- rev(seq( + round_any(range[1], step_size, f = ceiling), + round_any(est, step_size, f = floor), step_size + )) } else { - low <- rev(seq(round_any(range[1], step_size, f = ceiling), - round_any(est, step_size, f = floor) - step_size, step_size)) + low <- rev(seq( + round_any(range[1], step_size, f = ceiling), + round_any(est, step_size, f = floor) - step_size, step_size + )) } - if (est != round_any(est, step_size, f = ceiling)) { - high <- c(est, seq(round_any(est, step_size, f = ceiling), range[2], step_size)) + if (est != round_any(est, step_size, f = ceiling)) { + high <- c(est, seq(round_any(est, step_size, f = ceiling), range[2], step_size)) } else { - high <- c(seq(round_any(est, step_size, f = ceiling), range[2], step_size)) + high <- c(seq(round_any(est, step_size, f = ceiling), range[2], step_size)) } - vec <- c(low, high) - num <- sort(vec, index.return = TRUE)$ix - - profile <- r4ss::profile( - dir = profile_dir, - oldctlfile = "control.ss_new", - newctlfile = model_settings$newctlfile, - linenum = model_settings$linenum, - string = para, - profilevec = vec, - usepar = model_settings$usepar, - globalpar = model_settings$globalpar, - parlinenum = model_settings$parlinenum, - parstring = model_settings$parstring, - saveoutput = model_settings$saveoutput, - overwrite = model_settings$overwrite, - whichruns = model_settings$whichruns, - prior_check = model_settings$prior_check, - read_like = model_settings$read_like, - exe = model_settings$exe, - verbose = model_settings$verbose, - extras = model_settings$extras) - - # Save the output and the summary - name <- paste0("profile_", para) - vec_unordered <- vec - vec <- vec[num] - - profilemodels <- r4ss::SSgetoutput(dirvec = profile_dir, keyvec = num) - profilesummary <- r4ss::SSsummarize(biglist = profilemodels) - - profile_output <- list() - profile_output$mydir <- profile_dir - profile_output$para <- para - profile_output$name <- paste0("profile_", para) - profile_output$vec <- vec[num] - profile_output$model_settings <- model_settings - profile_output$profilemodels <- profilemodels - profile_output$profilesummary <- profilesummary - profile_output$rep <- rep - profile_output$vec_unordered <- vec - profile_output$num <- num - - save( - profile_dir, - para, - name, - vec, - vec_unordered, - model_settings, - profilemodels, - profilesummary, - rep, - num, - file = file.path(profile_dir, paste0(para, "_profile_output.Rdat")) - ) - - nwfscDiag::get_summary( - mydir = profile_dir, - name = paste0("profile_", para), - para = para, - # vec = vec[num], - vec = profilesummary$pars %>% - dplyr::filter(Label == para) %>% - dplyr::select(dplyr::starts_with("rep")) %>% - as.vector, - profilemodels = profilemodels, - profilesummary = profilesummary - ) - - nwfscDiag::profile_plot( - mydir = profile_dir, - para = para, - rep = rep, - profilesummary = profilesummary - ) - - message("Finished profile of ", para, ".") - - } #end parameter loop - -} \ No newline at end of file + vec <- c(low, high) + num <- sort(vec, index.return = TRUE)$ix + + profile <- r4ss::profile( + dir = profile_dir, + oldctlfile = "control.ss_new", + newctlfile = model_settings$newctlfile, + linenum = model_settings$linenum, + string = para, + profilevec = vec, + usepar = model_settings$usepar, + globalpar = model_settings$globalpar, + parlinenum = model_settings$parlinenum, + parstring = model_settings$parstring, + saveoutput = model_settings$saveoutput, + overwrite = model_settings$overwrite, + whichruns = model_settings$whichruns, + prior_check = model_settings$prior_check, + read_like = model_settings$read_like, + exe = model_settings$exe, + verbose = model_settings$verbose, + extras = model_settings$extras + ) + + # Save the output and the summary + name <- paste0("profile_", para) + vec_unordered <- vec + vec <- vec[num] + + profilemodels <- r4ss::SSgetoutput(dirvec = profile_dir, keyvec = num) + profilesummary <- r4ss::SSsummarize(biglist = profilemodels) + + profile_output <- list() + profile_output$mydir <- profile_dir + profile_output$para <- para + profile_output$name <- paste0("profile_", para) + profile_output$vec <- vec[num] + profile_output$model_settings <- model_settings + profile_output$profilemodels <- profilemodels + profile_output$profilesummary <- profilesummary + profile_output$rep <- rep + profile_output$vec_unordered <- vec + profile_output$num <- num + + save( + profile_dir, + para, + name, + vec, + vec_unordered, + model_settings, + profilemodels, + profilesummary, + rep, + num, + file = file.path(profile_dir, paste0(para, "_profile_output.Rdat")) + ) + + nwfscDiag::get_summary( + mydir = profile_dir, + name = paste0("profile_", para), + para = para, + # vec = vec[num], + vec = profilesummary$pars %>% + dplyr::filter(Label == para) %>% + dplyr::select(dplyr::starts_with("rep")) %>% + as.vector(), + profilemodels = profilemodels, + profilesummary = profilesummary + ) + + nwfscDiag::profile_plot( + mydir = profile_dir, + para = para, + rep = rep, + profilesummary = profilesummary + ) + + message("Finished profile of ", para, ".") + } # end parameter loop +} diff --git a/R/rerun_profile_vals.R b/R/rerun_profile_vals.R index d527a34..3225a80 100644 --- a/R/rerun_profile_vals.R +++ b/R/rerun_profile_vals.R @@ -1,174 +1,182 @@ #' Generate likelihood profiles #' To be called from the run_diagnostics function after creating #' the model settings using the get_settings function. -#' +#' #' #' @template mydir -#' @param para_name SS parameter name that the profile was run across. This is used to +#' @param para_name SS parameter name that the profile was run across. This is used to #' located the correct folder combined with the mydir function input (e.g. paste0(mydir, "_profile_", para_name)) #' @param run_num A single or vector of run numbers that you would like to rerun for convergence. #' This input needs to match the run number for the original profile (e.g., Report6.sso) that #' you would like to rerun. -#' @param data_file_nm SS data file name +#' @param data_file_nm SS data file name #' #' @author Chantel Wetzel. #' @export rerun_profile_vals <- function(mydir, - para_name, - run_num, - data_file_nm) -{ - - if(missing(mydir)){ - stop("Stop: Need to specify mydir.") - } - - if(missing(run_num)){ - stop("Stop: Need to specify run_num.") - } - - if(missing(para_name)){ - stop("Stop: Need to specify parameter name via parameter function input.") - } - para = para_name - - profile_dir = paste0(mydir, "_profile_", para_name) - temp_dir = file.path(profile_dir, "temp") - dir.create(temp_dir, showWarnings = FALSE) - - file.copy(file.path(profile_dir, 'ss.exe'), temp_dir, overwrite = TRUE) - file.copy(file.path(profile_dir, "control_modified.ss"), temp_dir) - file.copy(file.path(profile_dir, "control.ss_new"), temp_dir) - file.copy(file.path(profile_dir, data_file_nm), temp_dir) - file.copy(file.path(profile_dir, "starter.ss_new"), temp_dir) - file.copy(file.path(profile_dir, "forecast.ss_new"), temp_dir) - file.rename(file.path(temp_dir, "starter.ss_new"), file.path(temp_dir, "starter.ss")) - file.rename(file.path(temp_dir, "forecast.ss_new"), file.path(temp_dir, "forecast.ss")) - - # Use the SS_parlines funtion to ensure that the input parameter can be found - check_para <- r4ss::SS_parlines(dir = temp_dir, - verbose = FALSE, - active = FALSE)$Label == para - if( sum(check_para) == 0) { - stop(paste0( "The input profile_custom does not match a parameter in the control.ss_new file.")) - } - - load(file.path(profile_dir, paste0(para_name, "_profile_output.Rdat"))) - vec <- vec_unordered - like_check <- profilesummary$likelihoods[1,] - - # Change the control file name in the starter file - starter <- r4ss::SS_readstarter(file.path(temp_dir, 'starter.ss')) - starter$jitter_fraction <- 0.01 - starter$init_values_src <- model_settings$profile_init_values_src - # make sure the prior likelihood is calculated for non-estimated quantities - starter$prior_like <- model_settings$prior_like - r4ss::SS_writestarter(starter, dir = temp_dir, overwrite=TRUE) - - - for(i in run_num) { - setwd(temp_dir) - r4ss::SS_changepars( - ctlfile = "control_modified.ss", - newctlfile = "control_modified.ss", - strings = para, - newvals = vec[i], - estimate = FALSE - ) - system("ss -nohess") - - mod = r4ss::SS_output(dir = temp_dir, covar = FALSE, printstats = FALSE, verbose = FALSE) - like = mod$likelihoods_used[1, 1] - - # See if likelihood is lower than the original - and rerun if not - add <- 0.01 - if(like > like_check[i]) { - for(ii in 1:5) { - starter <- r4ss::SS_readstarter(file = file.path(temp_dir, 'starter.ss')) - starter$jitter_fraction <- add + starter$jitter_fraction - r4ss::SS_writestarter(starter, dir = temp_dir, overwrite = TRUE) - system("ss -nohess") - mod <- r4ss::SS_output(dir = temp_dir, covar = FALSE, printstats = FALSE, verbose = FALSE) - like <- mod$likelihoods_used[1, 1] - if( like < like_check[i] ){ break() } - } - } - - files = c("CompReport", "covar", "Report", "warning") - for(j in 1:length(files)){ - file.rename(paste0(files[j], ".sso"), - paste0(files[j], i, ".sso")) - file.copy(paste0(files[j], i, ".sso"), - profile_dir, overwrite = TRUE) - } - file.rename("ss.par", paste0("ss.par_", i, ".sso")) - file.copy(paste0("ss.par_", i, ".sso"), profile_dir, overwrite = TRUE) - - } - - profilemodels <- r4ss::SSgetoutput(dirvec = profile_dir, keyvec = num) - profilesummary <- r4ss::SSsummarize(biglist = profilemodels) - - protect = file.path(profile_dir, "protect") - dir.create(protect, showWarnings = FALSE) - save_files = c(paste0(para, "_trajectories_compare1_spawnbio.png"), - paste0(para, "_trajectories_compare3_Bratio.png"), - paste0("parameter_panel_", para, ".png"), - paste0("piner_panel_", para, ".png"), - paste0("profile_", para, "_parsonbounds.csv"), - paste0("profile_", para, "_quant_table.csv"), - paste0("profile_", para, "_results.csv"), - paste0(para, "_profile_output.Rdat")) - - for(i in 1:length(save_files)) { - file.copy(file.path(profile_dir, save_files[i]), - file.path(profile_dir, "protect")) - } - - name <- paste0("profile_", para) - vec_unordered <- vec - vec <- vec[num] - - profile_output <- list() - profile_output$mydir <- profile_dir - profile_output$para <- para - profile_output$name <- paste0("profile_", para) - profile_output$vec <- vec[num] - profile_output$model_settings <- model_settings - profile_output$profilemodels <- profilemodels - profile_output$profilesummary <- profilesummary - profile_output$rep <- rep - profile_output$vec_unordered <- vec - profile_output$num <- num - - save( - profile_dir, - para, - name, - vec, - vec_unordered, - model_settings, - profilemodels, - profilesummary, - rep, - num, - file = file.path(profile_dir, paste0(para, "_profile_output.Rdat")) - ) - - get_summary( - mydir = profile_dir, - name = paste0("profile_", para), - para = para, - vec = vec[num], - profilemodels = profilemodels, - profilesummary = profilesummary - ) - - profile_plot( - mydir = profile_dir, - para = para, - rep = rep, - profilesummary = profilesummary - ) - -} \ No newline at end of file + para_name, + run_num, + data_file_nm) { + if (missing(mydir)) { + stop("Stop: Need to specify mydir.") + } + + if (missing(run_num)) { + stop("Stop: Need to specify run_num.") + } + + if (missing(para_name)) { + stop("Stop: Need to specify parameter name via parameter function input.") + } + para <- para_name + + profile_dir <- paste0(mydir, "_profile_", para_name) + temp_dir <- file.path(profile_dir, "temp") + dir.create(temp_dir, showWarnings = FALSE) + + file.copy(file.path(profile_dir, "ss.exe"), temp_dir, overwrite = TRUE) + file.copy(file.path(profile_dir, "control_modified.ss"), temp_dir) + file.copy(file.path(profile_dir, "control.ss_new"), temp_dir) + file.copy(file.path(profile_dir, data_file_nm), temp_dir) + file.copy(file.path(profile_dir, "starter.ss_new"), temp_dir) + file.copy(file.path(profile_dir, "forecast.ss_new"), temp_dir) + file.rename(file.path(temp_dir, "starter.ss_new"), file.path(temp_dir, "starter.ss")) + file.rename(file.path(temp_dir, "forecast.ss_new"), file.path(temp_dir, "forecast.ss")) + + # Use the SS_parlines funtion to ensure that the input parameter can be found + check_para <- r4ss::SS_parlines( + dir = temp_dir, + verbose = FALSE, + active = FALSE + )$Label == para + if (sum(check_para) == 0) { + stop(paste0("The input profile_custom does not match a parameter in the control.ss_new file.")) + } + + load(file.path(profile_dir, paste0(para_name, "_profile_output.Rdat"))) + vec <- vec_unordered + like_check <- profilesummary$likelihoods[1, ] + + # Change the control file name in the starter file + starter <- r4ss::SS_readstarter(file.path(temp_dir, "starter.ss")) + starter$jitter_fraction <- 0.01 + starter$init_values_src <- model_settings$profile_init_values_src + # make sure the prior likelihood is calculated for non-estimated quantities + starter$prior_like <- model_settings$prior_like + r4ss::SS_writestarter(starter, dir = temp_dir, overwrite = TRUE) + + + for (i in run_num) { + setwd(temp_dir) + r4ss::SS_changepars( + ctlfile = "control_modified.ss", + newctlfile = "control_modified.ss", + strings = para, + newvals = vec[i], + estimate = FALSE + ) + system("ss -nohess") + + mod <- r4ss::SS_output(dir = temp_dir, covar = FALSE, printstats = FALSE, verbose = FALSE) + like <- mod$likelihoods_used[1, 1] + + # See if likelihood is lower than the original - and rerun if not + add <- 0.01 + if (like > like_check[i]) { + for (ii in 1:5) { + starter <- r4ss::SS_readstarter(file = file.path(temp_dir, "starter.ss")) + starter$jitter_fraction <- add + starter$jitter_fraction + r4ss::SS_writestarter(starter, dir = temp_dir, overwrite = TRUE) + system("ss -nohess") + mod <- r4ss::SS_output(dir = temp_dir, covar = FALSE, printstats = FALSE, verbose = FALSE) + like <- mod$likelihoods_used[1, 1] + if (like < like_check[i]) { + break() + } + } + } + + files <- c("CompReport", "covar", "Report", "warning") + for (j in 1:length(files)) { + file.rename( + paste0(files[j], ".sso"), + paste0(files[j], i, ".sso") + ) + file.copy(paste0(files[j], i, ".sso"), + profile_dir, + overwrite = TRUE + ) + } + file.rename("ss.par", paste0("ss.par_", i, ".sso")) + file.copy(paste0("ss.par_", i, ".sso"), profile_dir, overwrite = TRUE) + } + + profilemodels <- r4ss::SSgetoutput(dirvec = profile_dir, keyvec = num) + profilesummary <- r4ss::SSsummarize(biglist = profilemodels) + + protect <- file.path(profile_dir, "protect") + dir.create(protect, showWarnings = FALSE) + save_files <- c( + paste0(para, "_trajectories_compare1_spawnbio.png"), + paste0(para, "_trajectories_compare3_Bratio.png"), + paste0("parameter_panel_", para, ".png"), + paste0("piner_panel_", para, ".png"), + paste0("profile_", para, "_parsonbounds.csv"), + paste0("profile_", para, "_quant_table.csv"), + paste0("profile_", para, "_results.csv"), + paste0(para, "_profile_output.Rdat") + ) + + for (i in 1:length(save_files)) { + file.copy( + file.path(profile_dir, save_files[i]), + file.path(profile_dir, "protect") + ) + } + + name <- paste0("profile_", para) + vec_unordered <- vec + vec <- vec[num] + + profile_output <- list() + profile_output$mydir <- profile_dir + profile_output$para <- para + profile_output$name <- paste0("profile_", para) + profile_output$vec <- vec[num] + profile_output$model_settings <- model_settings + profile_output$profilemodels <- profilemodels + profile_output$profilesummary <- profilesummary + profile_output$rep <- rep + profile_output$vec_unordered <- vec + profile_output$num <- num + + save( + profile_dir, + para, + name, + vec, + vec_unordered, + model_settings, + profilemodels, + profilesummary, + rep, + num, + file = file.path(profile_dir, paste0(para, "_profile_output.Rdat")) + ) + + get_summary( + mydir = profile_dir, + name = paste0("profile_", para), + para = para, + vec = vec[num], + profilemodels = profilemodels, + profilesummary = profilesummary + ) + + profile_plot( + mydir = profile_dir, + para = para, + rep = rep, + profilesummary = profilesummary + ) +} diff --git a/R/run_diagnostics.R b/R/run_diagnostics.R index 5bb23d2..f05452e 100644 --- a/R/run_diagnostics.R +++ b/R/run_diagnostics.R @@ -2,16 +2,16 @@ #' 1) Jitter #' 2) Profiles across female m, steepness, and R0 #' 3) Retrospectives -#' +#' #' #' @template mydir #' @template model_settings -#' +#' #' @author Chantel Wetzel #' @return A vector of likelihoods for each jitter iteration. #' @export -run_diagnostics <- function(mydir, model_settings ){ +run_diagnostics <- function(mydir, model_settings) { # Check for Report file model_dir <- file.path(mydir, paste0(model_settings$base_name)) @@ -20,23 +20,23 @@ run_diagnostics <- function(mydir, model_settings ){ orig_dir <- getwd() setwd(model_dir) cat("Running model in directory:", getwd(), "\n") - run(dir = model_dir, - exe = model_settings$exe, - extras = model_settings$extras) + run( + dir = model_dir, + exe = model_settings$exe, + extras = model_settings$extras + ) setwd(orig_dir) } if ("retro" %in% model_settings$run) { retro_wrapper(mydir = mydir, model_settings = model_settings) - - } + } if ("profile" %in% model_settings$run) { - profile_wrapper(mydir = mydir, model_settings = model_settings) + profile_wrapper(mydir = mydir, model_settings = model_settings) } if ("jitter" %in% model_settings$run) { - jitter_wrapper(mydir = mydir, model_settings = model_settings) - } - + jitter_wrapper(mydir = mydir, model_settings = model_settings) + } }