This vignette includes an example of how to use parallel computing to conduct management strategy evaluation. Note: all code shown here is for demonstration purposes only; the results are for illustrative use and should not be interpreted as reflecting real-world stock assessments.
library(wham)
library(whamMSE)
main.dir = here::here()year_start <- 1 # starting year in the burn-in period
year_end <- 20 # end year in the burn-in period
MSE_years <- 3 # number of years in the feedback loop
# Note: no need to include MSE_years in simulation-estimation
info <- generate_basic_info(n_stocks = 2,
n_regions = 2,
n_indices = 2,
n_fleets = 2,
n_seasons = 4,
base.years = year_start:year_end,
n_feedback_years = MSE_years)
basic_info = info$basic_info # collect basic information
catch_info = info$catch_info # collect fleet catch information
index_info = info$index_info # collect survey information
F_info = info$F # collect fishing information
# see more details using ?generate_basic_infobasic_info <- generate_NAA_where(basic_info = basic_info, move.type = 2) # "bidirectional" movement
move <- generate_move(basic_info = basic_info, move.type = 2, move.rate = c(0.3,0.1), move.re = "constant")n_stocks <- as.integer(basic_info['n_stocks'])
n_regions <- as.integer(basic_info['n_regions'])
n_fleets <- as.integer(basic_info['n_fleets'])
n_indices <- as.integer(basic_info['n_indices'])
n_ages <- as.integer(basic_info['n_ages'])
# Selectivity Configuration
fleet_pars <- c(5,1)
index_pars <- c(2,1)
sel <- list(model=rep("logistic",n_fleets+n_indices),
initial_pars=c(rep(list(fleet_pars),n_fleets),rep(list(index_pars),n_indices)))M_em <- list(model="constant", initial_means=array(0.2, dim = c(n_stocks, n_regions, n_ages)))sigma <- "rec+1"
re_cor <- "iid"
ini.opt <- "equilibrium" # option <- c("age-specific-fe", "equilibrium")
Rec_sig <- 0.5 # (sigma for recruitment)
NAA_sig <- 0.5 # (sigma for NAA)
sigma_vals = array(NAA_sig, dim = c(n_stocks, n_regions, n_ages)) # n_stocks x n_regions x n_ages
sigma_vals[,,1] = Rec_sig
# Set initial NAA for each stock
log_N1 <- rep(10,n_stocks)
log_N1[1] <- log(exp(10)*2) # N1_stock1 is 2 times higher than N1_stock2
N1_pars <- generate_ini_N1(basic_info,ini.opt,log_N1)
# Set mean recruitment para. for each stock
mean_rec_par <- list()
for (i in 1:n_stocks) mean_rec_par[[i]] <- exp(log_N1[i])
NAA_re <- list(N1_model=rep(ini.opt,n_stocks),
sigma=rep(sigma,n_stocks),
cor=rep(re_cor,n_stocks),
recruit_model = 2, # rec random around the mean
recruit_pars = mean_rec_par,
sigma_vals = sigma_vals,
N1_pars = N1_pars)wham inputHere we use prepare_wham_input() function to generate a
wham input using the basic information we set above:
input <- prepare_wham_input(basic_info = basic_info,
selectivity = sel,
M = M,
NAA_re = NAA_re,
move = move,
catch_info = catch_info,
index_info = index_info,
F = F_info)random = input$random # check what processes are random effects
input$random = NULL # so inner optimization won't change simulated RE
om <- fit_wham(input, do.fit = F, do.brps = T, MakeADFun.silent = TRUE)
# Note: do.fit must be FALSE (no modeling fitting yet)assess.interval <- 3 # Note: assessment interval is 3 years, given the feedback period is 3 years, there will be only 1 assessment
base.years <- year_start:year_end # Burn-in period
first.year <- head(base.years,1)
terminal.year <- tail(base.years,1)
assess.years <- seq(terminal.year, tail(om$years,1)-assess.interval,by = assess.interval)
mods <- list() # Create a list to save MSE outputssub.dir = "Results"
dir.create(file.path(getwd(), sub.dir), recursive = TRUE)
library(doParallel)
library(foreach)
detectCores() # check how many cores availablecluster <- makeCluster(5)
registerDoParallel(cluster)
foreach (i = 1:5) %dopar% {
tryCatch({
library(wham)
library(whamMSE)
om_with_data <- update_om_fn(om, seed = 123+i, random = random)
n_stocks = n_regions = n_fleets = n_indices = 1
sel_em <- list(model=rep("logistic", n_fleets + n_indices),
initial_pars=c(rep(list(fleet_pars), n_fleets), rep(list(index_pars), n_indices)))
NAA_re_em <- list(N1_model="equilibrium", sigma="rec+1", cor="iid")
M_em <- list(model="constant", initial_means=array(0.2, dim = c(n_stocks, n_regions, n_ages)))
mod <- loop_through_fn(om = om_with_data,
em_info = info,
random = random,
M_em = M_em,
sel_em = sel_em,
NAA_re_em = NAA_re_em,
move_em = NULL,
em.opt = list(separate.em = TRUE, separate.em.type = 1, do.move = FALSE, est.move = FALSE),
aggregate_catch_info = list(n_fleets = 1, fleet_pointer = c(1,1),
use_catch_weighted_waa=TRUE, catch_Neff = 100, catch_cv = 0.1),
aggregate_index_info = list(n_indices = 1, index_pointer = c(1,1),
use_catch_weighted_waa=TRUE, index_Neff = 100, index_cv = 0.1),
assess_years = assess.years,
assess_interval = assess.interval,
base_years = base.years,
year.use = 20,
seed = 123+i)
saveRDS(mod, file.path(sub.dir, sprintf("Mod1_%03d.RDS", i)))
}, error = function(e) {
cat(sprintf("Error in iteration %d: %s\n", i, e$message))
})
}
stopCluster(cluster)cluster <- makeCluster(5)
registerDoParallel(cluster)
foreach (i = 1:5) %dopar% {
tryCatch({
library(wham)
library(whamMSE)
om_with_data <- update_om_fn(om, seed = 123+i, random = random)
n_stocks = n_regions = 1
n_fleets = n_indices = 2
sel_em <- list(model=rep("logistic", n_fleets + n_indices),
initial_pars=c(rep(list(fleet_pars), n_fleets), rep(list(index_pars), n_indices)))
NAA_re_em <- list(N1_model="equilibrium", sigma="rec+1", cor="iid")
M_em <- list(model="constant", initial_means=array(0.2, dim = c(n_stocks, n_regions, n_ages)))
mod <- loop_through_fn(om = om_with_data,
em_info = info,
random = random,
M_em = M_em,
sel_em = sel_em,
NAA_re_em = NAA_re_em,
move_em = NULL,
em.opt = list(separate.em = TRUE, separate.em.type = 2, do.move = FALSE, est.move = FALSE),
assess_years = assess.years,
assess_interval = assess.interval,
base_years = base.years,
year.use = 20,
seed = 123+i)
saveRDS(mod, file.path(sub.dir, sprintf("Mod2_%03d.RDS", i)))
}, error = function(e) {
cat(sprintf("Error in iteration %d: %s\n", i, e$message))
})
}
stopCluster(cluster)cluster <- makeCluster(5)
registerDoParallel(cluster)
foreach (i = 1:5) %dopar% {
tryCatch({
library(wham)
library(whamMSE)
om_with_data <- update_om_fn(om, seed = 123+i, random = random)
n_stocks = n_regions = n_fleets = n_indices = 1
sel_em <- list(model=rep("logistic", n_fleets + n_indices),
initial_pars=c(rep(list(fleet_pars), n_fleets), rep(list(index_pars), n_indices)))
NAA_re_em <- list(N1_model="equilibrium", sigma="rec+1", cor="iid")
M_em <- list(model="constant", initial_means=array(0.2, dim = c(n_stocks, n_regions, n_ages)))
mod <- loop_through_fn(om = om_with_data,
em_info = info,
random = random,
M_em = M_em,
sel_em = sel_em,
NAA_re_em = NAA_re_em,
move_em = NULL,
em.opt = list(separate.em = TRUE, separate.em.type = 3, do.move = FALSE, est.move = FALSE),
assess_years = assess.years,
assess_interval = assess.interval,
base_years = base.years,
year.use = 20,
seed = 123+i)
saveRDS(mod, file.path(sub.dir, sprintf("Mod3_%03d.RDS", i)))
}, error = function(e) {
cat(sprintf("Error in iteration %d: %s\n", i, e$message))
})
}
stopCluster(cluster)cluster <- makeCluster(5)
registerDoParallel(cluster)
foreach (i = 1:5) %dopar% {
tryCatch({
library(wham)
library(whamMSE)
om_with_data <- update_om_fn(om, seed = 123+i, random = random)
n_stocks = n_regions = n_fleets = n_indices = 2
sel_em <- list(model=rep("logistic", n_fleets + n_indices),
initial_pars=c(rep(list(fleet_pars), n_fleets), rep(list(index_pars), n_indices)))
NAA_re_em <- list(N1_model=rep("equilibrium", n_stocks),
sigma=rep("rec+1", n_stocks),
cor=rep("iid", n_stocks),
recruit_model = 2)
M_em <- list(model="constant", initial_means=array(0.2, dim = c(n_stocks, n_regions, n_ages)))
mod <- loop_through_fn(om = om_with_data,
em_info = info,
random = random,
M_em = M_em,
sel_em = sel_em,
NAA_re_em = NAA_re_em,
move_em = NULL,
age_comp_em = "multinomial",
em.opt = list(separate.em = FALSE, separate.em.type = NULL, do.move = FALSE, est.move = FALSE),
assess_years = assess.years,
assess_interval = assess.interval,
base_years = base.years,
year.use = 20,
seed = 123+i)
saveRDS(mod, file.path(sub.dir, sprintf("Mod4_%03d.RDS", i)))
}, error = function(e) {
cat(sprintf("Error in iteration %d: %s\n", i, e$message))
})
}
stopCluster(cluster)cluster <- makeCluster(5)
registerDoParallel(cluster)
foreach (i = 1:5) %dopar% {
tryCatch({
library(wham)
library(whamMSE)
om_with_data <- update_om_fn(om, seed = 123+i, random = random)
n_stocks = n_regions = n_fleets = n_indices = 2
sel_em <- list(model=rep("logistic", n_fleets + n_indices),
initial_pars=c(rep(list(fleet_pars), n_fleets), rep(list(index_pars), n_indices)))
NAA_re_em <- list(N1_model=rep("equilibrium", n_stocks),
sigma=rep("rec+1", n_stocks),
cor=rep("iid", n_stocks),
NAA_where = om$input$data$NAA_where)
M_em <- list(model="constant", initial_means=array(0.2, dim = c(n_stocks, n_regions, n_ages)))
mod <- loop_through_fn(om = om_with_data,
em_info = info,
random = random,
M_em = M_em,
sel_em = sel_em,
NAA_re_em = NAA_re_em,
move_em = move,
em.opt = list(separate.em = FALSE, separate.em.type = NULL, do.move = TRUE, est.move = FALSE),
assess_years = assess.years,
assess_interval = assess.interval,
base_years = base.years,
year.use = 20,
seed = 123+i)
saveRDS(mod, file.path(sub.dir, sprintf("Mod5_%03d.RDS", i)))
}, error = function(e) {
cat(sprintf("Error in iteration %d: %s\n", i, e$message))
})
}
stopCluster(cluster)