From e6bcb1a885949be004673dd3d8c44f345372ab81 Mon Sep 17 00:00:00 2001
From: Bai-Li-NOAA <bai.li@noaa.gov>
Date: Wed, 14 Aug 2024 12:30:59 -0400
Subject: [PATCH] add comparison of recruitment functions

---
 content/recruitment_functions_comparison.R | 161 +++++++++++++++++++
 content/run_pollock_tests_by_fleet.R       | 170 ---------------------
 2 files changed, 161 insertions(+), 170 deletions(-)
 create mode 100644 content/recruitment_functions_comparison.R
 delete mode 100644 content/run_pollock_tests_by_fleet.R

diff --git a/content/recruitment_functions_comparison.R b/content/recruitment_functions_comparison.R
new file mode 100644
index 0000000..5d964b7
--- /dev/null
+++ b/content/recruitment_functions_comparison.R
@@ -0,0 +1,161 @@
+source(file.path(getwd(), "content", "run_pollock_tests.R"))
+# Function to print arguments
+print_arguments <- function(...) {
+  cat("Following arguments were used:\n\n")
+  print(list(...))
+}
+
+# Option 1: Organize function parameters by individual recruitment model parameters --------
+setup_BevertonHolt_recruitment_opt1 <- function(log_rzero,
+                                                logit_steep,
+                                                log_sigma_recruit,
+                                                log_devs) {
+  recruitment_m <- methods::new(BevertonHoltRecruitment)
+
+  # Set up log_rzero parameters
+  recruitment_m$log_rzero$value <- log_rzero$value
+  recruitment_m$log_rzero$estimated <- log_rzero$estimated
+  recruitment_m$log_rzero$is_random_effect <- log_rzero$is_randome_effect
+
+  # Set up logit_steep parameters
+  recruitment_m$logit_steep$value <- logit_steep$value
+  recruitment_m$logit_steep$estimated <- logit_steep$estimated
+  recruitment_m$logit_steep$is_random_effect <- logit_steep$is_randome_effect
+
+  # Set up log_sigma_recruit parameters
+  recruitment_m$log_sigma_recruit$value <- log_sigma_recruit$value
+  recruitment_m$log_sigma_recruit$estimated <- log_sigma_recruit$estimated
+  recruitment_m$log_sigma_recruit$is_random_effect <- log_sigma_recruit$is_randome_effect
+
+  # Set up log_devs parameters
+  recruitment_m$log_devs <- log_devs$log_devs
+  recruitment_m$estimate_log_devs <- log_devs$estimate_log_devs
+
+  # Print the arguments used
+  print_arguments(
+    log_rzero = log_rzero,
+    logit_steep = logit_steep,
+    log_sigma_recruit = log_sigma_recruit,
+    log_devs = log_devs
+  )
+  
+  return(recruitment_m)
+}
+
+# Example: prepare inputs for option 1
+log_rzero <- list(
+  value = exp(parfinal$mean_log_recruit + log(1e9)),
+  estimated = TRUE,
+  is_randome_effect = FALSE
+)
+
+logit_steep <- list(
+  value = .99999,
+  estimated = FALSE,
+  is_randome_effect = FALSE
+)
+
+log_sigma_recruit <- logit_steep
+log_sigma_recruit$value <- parfinal$sigmaR
+
+log_devs <- list(
+  log_devs = parfinal$dev_log_recruit[-1],
+  estimate_log_devs = TRUE
+)
+
+# Apply option 1 function
+recruitment_opt1 <- setup_BevertonHolt_recruitment_opt1(
+  log_rzero,
+  logit_steep,
+  log_sigma_recruit,
+  log_devs
+)
+
+# Option 2: organize function parameters by fields of FIMS parameters (using a list of vectors) -----------------------------------------
+setup_BevertonHolt_recruitment_opt2 <- function(recruitment) {
+  r <- recruitment
+  recruitment_m <- methods::new(BevertonHoltRecruitment)
+  recruitment_m$log_sigma_recruit$value <- log(r$init_pars[3])
+  recruitment_m$log_rzero$value <- log(r$init_pars[1])
+  recruitment_m$log_rzero$is_random_effect <- FALSE
+  recruitment_m$log_rzero$estimated <- TRUE ## !(1 %in% r$fix_pars)
+  steep <- r$init_pars[2]
+  recruitment_m$logit_steep$value <- -log(1.0 - steep) + log(steep - 0.2)
+  recruitment_m$logit_steep$is_random_effect <- FALSE
+  recruitment_m$logit_steep$estimated <- FALSE ## !(2 %in% r$fix_pars)
+  recruitment_m$estimate_log_devs <- TRUE
+  recruitment_m$log_devs <- r$init_recdevs
+  print_arguments(recruitment = recruitment)
+  return(recruitment_m)
+}
+
+# Example: prepare inputs for option 2
+recruitment_inputs_opt2 <- list(
+  ## R0, steepness, sigmaR
+  init_pars = c(exp(parfinal$mean_log_recruit + log(1e9)), .99999, parfinal$sigmaR),
+  fix_pars = c(2:3),
+  init_recdevs = parfinal$dev_log_recruit[-1],
+  re = c("none", "none")
+)
+
+# Apply optoin 2 function
+recruitment_opt2 <- setup_BevertonHolt_recruitment_opt2(recruitment_inputs_opt2)
+
+# Option 3: Organize parameters by fields of FIMS parameters (using a list of named lists) -------------------------------------
+
+setup_BevertonHolt_recruitment_opt3 <- function(initial_values, is_estimated, is_random_effect) {
+  recruitment_m <- methods::new(BevertonHoltRecruitment)
+
+  # Set initial parameter values
+  recruitment_m$log_rzero$value <- initial_values$log_rzero
+  recruitment_m$logit_steep$value <- initial_values$logit_steep
+  recruitment_m$log_sigma_recruit$value <- initial_values$log_sigma_recruit
+  recruitment_m$log_devs <- initial_values$log_devs
+
+  # Set whether each parameter is estimated
+  recruitment_m$log_rzero$estimated <- is_estimated$log_rzero
+  recruitment_m$logit_steep$estimated <- is_estimated$logit_steep
+  recruitment_m$log_sigma_recruit$estimated <- is_estimated$log_sigma_recruit
+  recruitment_m$estimate_log_devs <- is_estimated$log_devs
+
+  # Set whether each parameter has random effects
+  recruitment_m$log_rzero$is_random_effect <- is_random_effect$log_rzero
+  recruitment_m$logit_steep$is_random_effect <- is_random_effect$logit_steep
+  recruitment_m$log_sigma_recruit$estimated <- is_estimated$log_sigma_recruit
+
+  print_arguments(
+    initial_values = initial_values,
+    is_estimated = is_estimated,
+    is_random_effect = is_random_effect
+  )
+  
+  return(recruitment_m)
+}
+
+# Example: prepare inputs for option 3
+recruitment_inputs_opt3 <- list(
+  initial_values = list(
+    log_rzero = parfinal$mean_log_recruit + log(1e9),
+    logit_steep = -log(1.0 - .99999) + log(.99999 - 0.2),
+    log_sigma_recruit = log(parfinal$sigmaR),
+    log_devs = parfinal$dev_log_recruit[-1]
+  ),
+  is_estimated = list(
+    log_rzero = TRUE,
+    logit_steep = FALSE,
+    log_sigma_recruit = FALSE,
+    log_devs = TRUE
+  ),
+  is_random_effect = list(
+    log_rzero = FALSE,
+    logit_steep = FALSE,
+    log_sigma_recruit = FALSE
+  )
+)
+
+# Apply option 3 function
+recruitment_opt3 <- setup_BevertonHolt_recruitment_opt3(
+  initial_values = recruitment_ctl$initial_values,
+  is_estimated = recruitment_ctl$is_estimated,
+  is_random_effect = recruitment_ctl$is_random_effect
+)
diff --git a/content/run_pollock_tests_by_fleet.R b/content/run_pollock_tests_by_fleet.R
deleted file mode 100644
index e479404..0000000
--- a/content/run_pollock_tests_by_fleet.R
+++ /dev/null
@@ -1,170 +0,0 @@
-## define the dimensions and global variables
-# library(FIMS)
-# 
-# source(file.path(getwd(), "content", "R", "pk_prepare_FIMS_inputs_by_fleet.R"))
-# 
-# years <- 1970:2023
-# nyears <- length(years)
-# nseasons <- 1
-# nages <- 10
-# ages <- 1:nages
-# source(file.path(getwd(), "content", "R", "pk_prepare_dat.R"))
-
-# test R functions: approach 1 --------------------------------------------
-clear()
-clear_logs()
-
-
-
-## make FIMS model
-success <- CreateTMBModel()
-parameters <- list(p = get_fixed())
-obj3 <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS", silent = TRUE)
-## report values for the two models
-rep3 <- obj3$report()
-
-
-clear()
-clear_logs()
-# source("R/pk_prepare_FIMS_inputs.R")
-estimate_fish_selex <- TRUE
-estimate_survey_selex <- TRUE
-estimate_q2 <- TRUE
-estimate_q3 <- TRUE
-estimate_q6 <- TRUE
-estimate_F <- TRUE
-estimate_recdevs <- TRUE
-source(file.path(getwd(), "content", "R", "pk_prepare_FIMS_inputs.R"))
-## make FIMS model
-success <- CreateTMBModel()
-parameters <- list(p = get_fixed())
-obj1 <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS", silent = TRUE)
-rep1 <- obj1$report()
-
-## Parameters are in a different order but seem to match
-opt1 <- with(obj1, nlminb(par, fn, gr))
-opt3 <- with(obj3, nlminb(par, fn, gr))
-clear()
-clear_logs()
-
-
-## Test that the two models are identical
-all.equal(rep1,rep3)
-all.equal(obj1$report(opt1$par), obj3$report(opt3$par))
-sum(obj1$par)-sum(obj3$par)
-
-print(cbind(sort(opt1$par)- sort(opt3$par)))
-
-# test R functions: approach 2 --------------------------------------------
-# selectivity_double_logistic <- set_selectivity(
-#   form = "DoubleLogisticSelectivity",
-#   initial_values = c(
-#     parfinal$inf1_fsh_mean, exp(parfinal$log_slp1_fsh_mean),
-#     parfinal$inf2_fsh_mean, exp(parfinal$log_slp2_fsh_mean)
-#   ),
-#   is_estimated = rep(TRUE, 4),
-#   is_random_effect = rep(FALSE, 4)
-# )
-# 
-# selectivity_logistic <- set_selectivity(
-#   form = "LogisticSelectivity",
-#   initial_values = c(
-#     parfinal$inf1_fsh_mean, exp(parfinal$log_slp1_fsh_mean)
-#   ),
-#   is_estimated = rep(TRUE, 2),
-#   is_random_effect = rep(FALSE, 2)
-# )
-# 
-# selectivity <- vector(mode = "list", length = 4)
-# selectivity[[1]] <- selectivity[[2]] <- selectivity_double_logistic
-# selectivity[[3]] <- selectivity_logistic
-# selectivity[[4]] <- selectivity_double_logistic # Order matters because the selectivity modules assign unique IDs.
-# 
-# fleet <- set_fleet(
-#   data = age_frame,
-#   is_estimated_obs_error = list(FALSE, FALSE, FALSE, FALSE),
-#   selectivity_ctl = NULL,
-#   selectivity_module_list = selectivity,
-#   Fmort_ctl = list(
-#     is_survey = list(FALSE, TRUE, TRUE, TRUE),
-#     estimate_F = list(TRUE, FALSE, FALSE, FALSE),
-#     random_F = list(FALSE, FALSE, FALSE, FALSE)
-#   ),
-#   catchability_ctl = list(
-#     log_q = list(
-#       0, parfinal$log_q2_mean,
-#       parfinal$log_q3_mean, parfinal$inf1_srv6
-#     ),
-#     estimate_q = list(FALSE, TRUE, TRUE, TRUE),
-#     random_q = list(FALSE, FALSE, FALSE, FALSE)
-#   )
-# )
-# 
-# # Population module
-# # recruitment
-# recruitment <- set_recruitment(
-#   form = "BevertonHoltRecruitment",
-#   initial_values = list(
-#     log_rzero = parfinal$mean_log_recruit + log(1e9),
-#     logit_steep = -log(1.0 - .99999) + log(.99999 - 0.2),
-#     log_sigma_recruit = log(parfinal$sigmaR),
-#     log_devs = parfinal$dev_log_recruit[-1]
-#   ),
-#   is_estimated = list(
-#     log_rzero = TRUE,
-#     logit_steep = FALSE,
-#     log_sigma_recruit = FALSE,
-#     log_devs = TRUE
-#   ),
-#   is_random_effect = list(
-#     log_rzero = FALSE,
-#     logit_steep = FALSE,
-#     log_sigma_recruit = FALSE
-#   )
-# )
-# 
-# ## growth  -- assumes single WAA vector for everything, based on
-# ## Srv1 above
-# growth <- set_growth(
-#   data = population_data,
-#   form = "EWAAgrowth",
-#   initial_values = list(weights = pkinput$dat$wt_srv1)
-# )
-# ## NOTE: FIMS assumes SSB calculated at the start of the year, so
-# ## need to adjust ASAP to do so as well for now, need to make
-# ## timing of SSB calculation part of FIMS later
-# ## maturity
-# ## NOTE: for now tricking FIMS into thinking age 0 is age 1, so need to adjust A50 for maturity because FIMS calculations use ages 0-5+ instead of 1-6
-# maturity <- set_maturity(
-#   form = maturity_ctl$form,
-#   initial_values = maturity_ctl$initial_values,
-#   is_estimated = maturity_ctl$is_estimated,
-#   is_random_effect = maturity_ctl$is_random_effect
-# )
-# 
-# # population
-# population <- set_population(
-#   data = population_data,
-#   initial_values = population_initial_values,
-#   is_estimated = population_is_estimated,
-#   maturity = maturity,
-#   growth = growth,
-#   recruitment = recruitment
-# )
-# 
-# ## make FIMS model
-# success <- CreateTMBModel()
-# parameters <- list(p = get_fixed())
-# obj4 <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS", silent = TRUE)
-# ## report values for the two models
-# rep4 <- obj4$report()
-# 
-# ## Parameters are in a different order but seem to match
-# opt4 <- with(obj4, nlminb(par, fn, gr))
-# clear()
-# clear_logs()
-# # rep4
-# # rm(list = ls())
-# 
-# ## Test that the two models are identical
-# all.equal(rep3, rep4)