added checks for FXSPR and F from catch in projections and vectorized…

… which_F_age

R/compare_wham_models.R

@@ -43,7 +43,7 @@
 #'     \item{\code{$alpha}}{scalar, (1-alpha)\% confidence intervals will be plotted. Default = 0.05 for 95\% CI.}
 #'     \item{\code{$ages.lab}}{vector, overwrite model age labels.}
 #'     \item{\code{$kobe.yr}}{integer, which year to use in Kobe plot (relative status). Default = terminal model year.}
-#'     \item{\code{$M.age}}{integer, which age to use in M time-series plot. Default = `data$which_F_age` (age of F to use for full total F).}
+#'     \item{\code{$M.age}}{integer, which age to use in M time-series plot. Default = `max(data$which_F_age)` (max age of F to use for full total F).}
 #'     \item{\code{$return.ggplot}}{T/F, return a list of ggplot2 objects for later modification? Default = TRUE.}
 #'   }
 #'
@@ -185,7 +185,7 @@ Returning AIC/rho table for WHAM models only.
       if(is.null(plot.opts$ages.lab)) plot.opts$ages.lab <- paste0(1:dim(x[[1]]$MAA)[2], c(rep("",dim(x[[1]]$MAA)[2]-1),"+"))
     }
     if(is.null(plot.opts$M.age)){
-      if(!no.wham) plot.opts$M.age <- wham.mods[[1]]$env$data$which_F_age
+      if(!no.wham) plot.opts$M.age <- max(wham.mods[[1]]$env$data$which_F_age)
       if(is.null(plot.opts$M.age)) plot.opts$M.age <- dim(x[[1]]$MAA)[2]
     }
 

R/fit_wham.R

@@ -100,6 +100,8 @@ fit_wham = function(input, n.newton = 3, do.sdrep = TRUE, do.retro = TRUE, n.pee
     btime <- Sys.time()
     mod <- fit_tmb(mod, n.newton = n.newton, do.sdrep = do.sdrep, do.check = do.check, save.sdrep = save.sdrep)
     mod$runtime = round(difftime(Sys.time(), btime, units = "mins"),2) # don't count retro or proj in runtime
+    mod = check_FXSPR(mod)
+    if(mod$env$data$do_proj==1) mod = check_projF(mod) #projections added.
 
     # retrospective analysis
     if(do.retro){
@@ -138,3 +140,76 @@ fit_wham = function(input, n.newton = 3, do.sdrep = TRUE, do.retro = TRUE, n.pee
 
   return(mod)
 }
+
+
+check_FXSPR = function(mod)
+{
+  percentSPR_out = exp(mod$rep$log_SPR_FXSPR - mod$rep$log_SPR0)
+  ind = which(round(percentSPR_out,4) != round(mod$env$data$percentSPR/100,4))
+  if(length(ind))
+  {
+    for(i in 1:2)
+    {
+      if(mod$env$data$n_years_proj) years = mod$years_full
+      else years = mod$years
+      print(years)
+      print(ind)
+      redo_SPR_years = years[ind]
+      print(redo_SPR_years)
+      print(paste(redo_SPR_years, collapse = ","))
+      warning(paste0("Changing initial values for estimating FXSPR for years ", paste(redo_SPR_years, collapse = ","), "."))
+      mod$env$data$FXSPR_init[ind] = mod$env$data$FXSPR_init[ind]*0.5
+      mod$fn(mod$opt$par)
+      mod$rep = mod$report(mod$opt$par)
+      percentSPR_out = exp(mod$rep$log_SPR_FXSPR - mod$rep$log_SPR0)
+      ind = which(round(percentSPR_out,4) != round(mod$env$data$percentSPR/100,4))
+      if(!length(ind)) break
+    }
+  }
+  if(length(ind)) warning(paste0("Still bad initial values and estimates of FXSPR for years ", paste(years[ind], collapse = ","), "."))
+  return(mod)     
+}
+
+check_projF = function(mod)
+{
+  proj_F_opt = mod$env$data$proj_F_opt
+  ind = which(proj_F_opt == 3) #FXSPR
+  if(length(ind))
+  {
+    y = mod$env$data$n_years_model + ind
+    bad = which(round(exp(mod$rep$log_FXSPR[y]),4) != round(mod$rep$FAA_tot[cbind(y,mod$env$data$which_F_age)],4))
+    if(length(bad))
+    {
+      redo_SPR_years = mod$years_full[y[bad]]
+      warning(paste0("Changing initial values for estimating FXSPR used to define F in projection years ", paste(redo_SPR_years, collapse = ","), "."))
+      mod$env$data$F_init_proj[ind[bad]] = mod$env$data$FXSPR_init_proj[y[bad]]
+      mod$fn(mod$opt$par)
+      mod$rep = mod$report(mod$opt$par)
+      bad = which(round(exp(mod$rep$log_FXSPR[y]),4) != round(mod$rep$FAA_tot[cbind(y,mod$env$data$which_F_age)],4))
+    }
+  }
+  y_bad_FXSPR = mod$years_full[y[bad]]
+  if(length(bad)) warning(paste0("Still bad initial values and estimates of FXSPR used to define F in projection years ", paste(y_bad_FXSPR, collapse = ","), "."))
+  ind = which(proj_F_opt == 5) #Find F from catch
+  if(length(ind))
+  {
+    y = mod$env$data$n_years_model + ind
+    bad = which(round(mod$env$data$data$proj_Fcatch[ind],4) != round(sum(mod$rep$pred_catch[y,]),4))
+    if(length(bad))
+    {
+      for(i in 1:2)
+      {
+        redo_Catch_years = mod$years_full[y[bad]]
+        warning(paste0("Changing initial values for finding F from Catch in projection years ", paste(redo_Catch_years, collapse = ","), , "."))
+        mod$env$data$F_init_proj[ind[bad]] = mod$env$data$F_init_proj[ind[bad]]*0.5
+        mod$fn(mod$opt$par)
+        mod$rep = mod$report(mod$opt$par)
+        bad = which(round(mod$env$data$data$proj_Fcatch[ind],4) != round(sum(mod$rep$pred_catch[y,]),4))
+        if(!length(bad)) break
+      }
+    }
+  }
+  y_bad_Fcatch = mod$years_full[y[bad]]
+  if(length(bad)) warning(paste0("Still bad initial values for finding F from Catch in projection years ", paste(y_bad_Fcatch, collapse = ","), , "."))
+  return(mod)
+}

R/prepare_projection.R

@@ -136,9 +136,14 @@ prepare_projection = function(model, proj.opts)
 #         capture.output(cat(c("use.last.F","use.avg.F","use.FXSPR","proj.F","proj.catch")[data$proj_F_opt])),"",sep='\n'))      
 #     }
 
+  data$FXSPR_init = c(data$FXSPR_init,rep(data$FXSPR_init[data$n_years_model], data$n_years_proj))
+  data$FMSY_init_proj = c(data$FMSY_init,rep(data$FMSY_init[data$n_years_model], data$n_years_proj))
+  data$F_proj_init = rep(0.1, data$n_years_proj)
+  data$F_proj_init[which(data$proj_F_opt == 3)] = data$FXSPR_init[which(data$proj_F_opt == 3)]
+  data$F_proj_init[which(data$proj_F_opt == 6)] = data$FMSY_init[which(data$proj_F_opt == 6)]
   #define age for full F in projections
   FAA_proj = colMeans(rbind(model$rep$FAA_tot[avg.yrs.ind,]))
-  data$which_F_age = which(FAA_proj == max(FAA_proj))[1]
+  data$which_F_age = c(data$which_F_age, rep(which(FAA_proj == max(FAA_proj))[1], data$n_years_proj))
 
   # modify data objects for projections (pad with average over avg.yrs): mature, fracyr_SSB, waa
   avg_cols = function(x) apply(x, 2, mean, na.rm=TRUE)

R/prepare_wham_input.R

@@ -539,7 +539,7 @@ without changing ASAP file, specify M$initial_means.")
 
   # data$recruit_model = 2 #random about mean
   data$N1_model = 0 #0: just age-specific numbers at age
-  data$which_F_age = data$n_ages #plus group by default used to define full F and F RP IN projections, only. prepare_projection changes it to properly define selectivity for projections.
+  data$which_F_age = rep(data$n_ages,data$n_years_model) #plus group by default used to define full F and F RP IN projections, only. prepare_projection changes it to properly define selectivity for projections.
   data$use_steepness = 0 #use regular SR parameterization by default, steepness still can be estimated as derived par.
   data$bias_correct_pe = 1 #bias correct log-normal process errors?
   data$bias_correct_oe = 1 #bias correct log-normal observation errors?
@@ -552,6 +552,7 @@ without changing ASAP file, specify M$initial_means.")
   # data$XSPR_R_opt = 3 #1(3): use annual R estimates(predictions) for annual SSB_XSPR, 2(4): use average R estimates(predictions). See next line for years to average over.
   data$XSPR_R_opt = 2 # default = use average R estimates
   data$XSPR_R_avg_yrs = 1:data$n_years_model #model year indices (TMB, starts @ 0) to use for averaging recruitment when defining SSB_XSPR (if XSPR_R_opt = 2,4)
+
   model_years <- asap3$year1 + 1:asap3$n_years - 1
 
   # --------------------------------------------------------------------------------
@@ -918,6 +919,10 @@ Ex: ",ecov$label[i]," in ",years[1]," affects ", c('recruitment','M')[data$Ecov_
   data$proj_M_opt <- 0
   data$logR_mean <- 0 # only used for SCAA projections
   data$logR_sd <- 0 # only used for SCAA projections
+  data$FXSPR_init = rep(0.1, data$n_years_model + data$n_years_proj)
+  data$FMSY_init = rep(0.1, data$n_years_model + data$n_years_proj)
+  data$F_proj_init = 0
+
   # data$obsvec[data$keep_I[data$use_indices==1]+1] - log(data$agg_indices[data$use_indices==1])
   # data$obsvec[data$keep_E[data$Ecov_use_obs==1]+1] - data$Ecov_obs[data$Ecov_use_obs==1]
 

R/prepare_wham_om_input.R

@@ -485,7 +485,7 @@ without changing ASAP file, specify M$initial_means.")
 
   # data$recruit_model = 2 #random about mean
   data$N1_model = 0 #0: just age-specific numbers at age
-  data$which_F_age = data$n_ages #plus group by default used to define full F and F RP IN projections, only. prepare_projection changes it to properly define selectivity for projections.
+  data$which_F_age = rep(data$n_ages,data$n_years_model) #plus group by default used to define full F and F RP IN projections, only. prepare_projection changes it to properly define selectivity for projections.
   data$use_steepness = basic_info$use_steepness #0: regular SR parameterization by default, steepness still can be estimated as derived par.
   data$bias_correct_pe = 0 #bias correct log-normal process errors?
   data$bias_correct_oe = 0 #bias correct log-normal observation errors?

R/prepare_wham_om_proj.R

@@ -123,7 +123,7 @@ prepare_wham_om_proj = function(om_input, proj.opts)
 
   #define age for full F in projections
   FAA_proj = colMeans(rbind(rep$FAA_tot[avg.yrs.ind,]))
-  data$which_F_age = which(FAA_proj == max(FAA_proj))[1]
+  data$which_F_age = c(data$which_F_age, rep(which(FAA_proj == max(FAA_proj))[1], data$n_years_proj))
 
   # modify data objects for projections (pad with average over avg.yrs): mature, fracyr_SSB, waa
   avg_cols = function(x) apply(x, 2, mean, na.rm=TRUE)

R/project_wham.R

@@ -125,3 +125,5 @@ project_wham = function(model, proj.opts=list(n.yrs=3, use.last.F=TRUE, use.avg.
   }
   return(mod)
 }
+
+

R/update_wham_input.R

@@ -353,7 +353,7 @@ Ex: ",ecov$label[i]," in ",years[1]," affects ", c('recruitment','M')[data$Ecov_
   #par$log_F1 = rep(-2, data$n_fleets)
   F = matrix(NA, nym,data$n_fleets)
   par$F_devs = matrix(0, nym-1, data$n_fleets)
-  for(i in 1:data$n_fleets) par$F_devs[,i] = diff(log(simres$FAA[1:nym,i,data$which_F_age]))
+  for(i in 1:data$n_fleets) par$F_devs[,i] = diff(log(simres$FAA[cbind(1:nym,i,data$which_F_age[1:nym])]))
 
   #par$F_devs = matrix(0, data$n_years_model-1, data$n_fleets)
   #if(data$N1_model == 1) par$log_N1_pars = c(10,log(0.1))

inst/example_scripts/ex8_compare_asap.R

@@ -4,6 +4,7 @@
 # devtools::install_github("timjmiller/wham", ref="devel", dependencies=TRUE)
 
 # devtools::load_all("~/Documents/wham/") # test locally (not yet pushed to github)
+# devtools::load_all("~/work/wham/wham")
 # devtools::install_github("timjmiller/wham", ref="devel", dependencies=TRUE)
 # devtools::install_github("timjmiller/wham", dependencies=TRUE)
 library(tidyverse)
@@ -15,6 +16,7 @@ main.dir <- "/home/bstock/Documents/wham/sandbox/wham_testing/runall-20210330"
 wham.dir <- find.package("wham")
 asap.dir <- file.path(wham.dir,"extdata","BASE_3")
 # asap.dir <- "/home/bstock/Documents/wham/inst/extdata/BASE_3"
+# asap.dir <- "~/work/wham/wham/inst/extdata/BASE_3"
 
 if(!dir.exists(main.dir)) dir.create(main.dir)
 write.dir <- file.path(main.dir,"ex8")