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.

1. Load “WHAM” and “whamMSE”:

library(wham)
library(whamMSE)

main.dir = here::here()

2. Generate Baisc Information

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_info

3. Specify Movement

basic_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")

4. Configure Selecitvity

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)))

5. Configure Natural Mortality

M_em <- list(model="constant", initial_means=array(0.2, dim = c(n_stocks, n_regions, n_ages)))

6. Configure numbers-at-age (NAA)

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)

7. Generate wham input

Here 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)

8. Generate operating model

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)

9. Specify Assessment Interval

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 outputs

10. Create A Folder to Save Results

sub.dir = "Results"
dir.create(file.path(getwd(), sub.dir), recursive = TRUE)

library(doParallel)
library(foreach)

detectCores() # check how many cores available

11. One-region Panmictic Assessment Model

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 = 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)

12. Spatially Implicit Assessment Model

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)

13. Separate Panmictic Assessment Models

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)

14. Spatially Disaggregated without Movement

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)

15. Spatially Explicit with Movement Fixed as Known

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)