From a03d2fe224fa1192ccb376e2cc99c8b2631fc7f4 Mon Sep 17 00:00:00 2001
From: Bai-Li-NOAA <bai.li@noaa.gov>
Date: Thu, 1 Aug 2024 09:49:35 -0400
Subject: [PATCH] add code to prepare FIMS inputs by fleet

---
 .../run-pollock-tests-by-process.yml          |  39 --
 ...sts-by-fleet.yml => run-pollock-tests.yml} |  19 +-
 content/R/pk_prepare_FIMS_inputs.R            | 495 ++++++++--------
 content/R/pk_prepare_FIMS_inputs_by_fleet.R   | 531 +++++++++++++-----
 content/R/pk_prepare_FIMS_inputs_by_process.R |  18 +-
 content/R/pk_prepare_dat.R                    |  19 -
 content/run_pollock_tests.R                   |  90 +--
 content/run_pollock_tests_by_fleet.R          | 210 ++-----
 8 files changed, 737 insertions(+), 684 deletions(-)
 delete mode 100644 .github/workflows/run-pollock-tests-by-process.yml
 rename .github/workflows/{run-pollock-tests-by-fleet.yml => run-pollock-tests.yml} (74%)

diff --git a/.github/workflows/run-pollock-tests-by-process.yml b/.github/workflows/run-pollock-tests-by-process.yml
deleted file mode 100644
index c5afddc..0000000
--- a/.github/workflows/run-pollock-tests-by-process.yml
+++ /dev/null
@@ -1,39 +0,0 @@
-on:
-  workflow_dispatch:
-
-name: Run pollock tests (by process)
-
-jobs:
-  build-deploy:
-    runs-on: macos-latest
-    env:
-          GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
-
-    steps:
-      - name: Check out repository
-        uses: actions/checkout@v2 
-        
-      - name: Set up R (needed for Rmd)
-        uses: r-lib/actions/setup-r@v2
-
-      - name: Run pollock tests 
-        run: |
-          # Names of required packages
-          packages <- c("TMB")
-
-          # Install packages not yet installed
-          installed_packages <- packages %in% rownames(installed.packages())
-          if (any(installed_packages == FALSE)) {
-            install.packages(packages[!installed_packages], repos = "http://cran.us.r-project.org")
-          }
-          
-          install.packages("FIMS", repos = c("https://noaa-fims.r-universe.dev", "https://cloud.r-project.org"))
-
-          # Load packages
-          invisible(lapply(packages, library, character.only = TRUE))
-
-          source(file.path(getwd(), "content", "run_pollock_tests.R"))
-        shell: Rscript {0}
-      
-      - name: Upload artifact
-        uses: actions/upload-pages-artifact@v3
diff --git a/.github/workflows/run-pollock-tests-by-fleet.yml b/.github/workflows/run-pollock-tests.yml
similarity index 74%
rename from .github/workflows/run-pollock-tests-by-fleet.yml
rename to .github/workflows/run-pollock-tests.yml
index abbc654..cf22bf8 100644
--- a/.github/workflows/run-pollock-tests-by-fleet.yml
+++ b/.github/workflows/run-pollock-tests.yml
@@ -1,7 +1,8 @@
+name: Run pollock tests
+
 on:
   workflow_dispatch:
-
-name: Run pollock tests (by fleet)
+  push:
 
 jobs:
   build-deploy:
@@ -19,7 +20,7 @@ jobs:
       - name: Run pollock tests 
         run: |
           # Names of required packages
-          packages <- c("dplyr", "tidyr", "ggplot2", "TMB", "reshape2", "here", "remotes", "lubridate")
+          packages <- c("dplyr", "tidyr", "ggplot2", "TMB", "reshape2", "here", "devtools", "lubridate")
 
           # Install packages not yet installed
           installed_packages <- packages %in% rownames(installed.packages())
@@ -27,13 +28,11 @@ jobs:
             install.packages(packages[!installed_packages], repos = "http://cran.us.r-project.org")
           }
           
-          install.packages("FIMS", repos = c("https://noaa-fims.r-universe.dev", "https://cloud.r-project.org"))
-
+          # install.packages("FIMS", repos = c("https://noaa-fims.r-universe.dev", "https://cloud.r-project.org"))
+          devtools::install_github("NOAA-FIMS/FIMS")
+          
           # Load packages
           invisible(lapply(packages, library, character.only = TRUE))
 
-          source(file.path(getwd(), "content", "run_pollock_tests_by_fleet.R"))
-        shell: Rscript {0}
-      
-      - name: Upload artifact
-        uses: actions/upload-pages-artifact@v3
+          source(file.path(getwd(), "content", "run_pollock_tests.R"))
+        shell: Rscript {0}
\ No newline at end of file
diff --git a/content/R/pk_prepare_FIMS_inputs.R b/content/R/pk_prepare_FIMS_inputs.R
index 02840e9..f8cce47 100644
--- a/content/R/pk_prepare_FIMS_inputs.R
+++ b/content/R/pk_prepare_FIMS_inputs.R
@@ -1,159 +1,164 @@
+# ## build a FIMS and PK data set that match
+# ##  need to fill missing years with -999 so it's ignored in FIMS
+# ind2 <- 0 * pkfit0$rep$Eindxsurv2 - 999
+# ind2[which(years %in% fimsdat$srvyrs2)] <- fimsdat$indxsurv2
+# CV2 <- rep(1, length = nyears) # actually SE in log space
+# CV2[which(years %in% fimsdat$srvyrs2)] <- fimsdat$indxsurv_log_sd2
+# paa2 <- pkfit0$rep$Esrvp2 * 0 - 999
+# paa2[which(years %in% fimsdat$srv_acyrs2), ] <- fimsdat$srvp2
+# Npaa2 <- rep(1, nyears)
+# Npaa2[which(years %in% fimsdat$srv_acyrs2)] <- fimsdat$multN_srv2
+# 
+# ind3 <- 0 * pkfit0$rep$Eindxsurv3 - 999
+# ind3[which(years %in% fimsdat$srvyrs3)] <- fimsdat$indxsurv3
+# CV3 <- rep(1, length = nyears) # actually SE in log space
+# CV3[which(years %in% fimsdat$srvyrs3)] <- fimsdat$indxsurv_log_sd3
+# paa3 <- pkfit0$rep$Esrvp3 * 0 - 999
+# paa3[which(years %in% fimsdat$srv_acyrs3), ] <- fimsdat$srvp3
+# Npaa3 <- rep(1, nyears)
+# Npaa3[which(years %in% fimsdat$srv_acyrs3)] <- fimsdat$multN_srv3
+# 
+# ind6 <- 0 * pkfit0$rep$Eindxsurv6 - 999
+# ind6[which(years %in% fimsdat$srvyrs6)] <- fimsdat$indxsurv6
+# CV6 <- rep(1, length = nyears) # actually SE in log space
+# CV6[which(years %in% fimsdat$srvyrs6)] <- fimsdat$indxsurv_log_sd6
+# paa6 <- pkfit0$rep$Esrvp6 * 0 - 999
+# paa6[which(years %in% fimsdat$srv_acyrs6), ] <- fimsdat$srvp6
+# Npaa6 <- rep(1, nyears)
+# Npaa6[which(years %in% fimsdat$srv_acyrs6)] <- fimsdat$multN_srv6
+# 
+# ## repeat with fish catch at age, using expected in missing years
+# caa <- pkfit0$rep$Ecatp * 0 - 999
+# caa[which(years %in% fimsdat$fshyrs), ] <- fimsdat$catp
+# Ncaa <- rep(1, nyears)
+# Ncaa[which(years %in% fimsdat$fshyrs)] <- fimsdat$multN_fsh
+# 
+# 
+# 
+# ## put into fims friendly form
+# res <- data.frame(
+#   type = character(),
+#   name = character(),
+#   age = integer(),
+#   datestart = character(),
+#   dateend = character(),
+#   value = double(),
+#   unit = character(),
+#   uncertainty = double()
+# )
+# landings <- data.frame(
+#   type = "landings",
+#   name = "fleet1",
+#   age = NA,
+#   datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"),
+#   dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"),
+#   value = as.numeric(fimsdat$cattot) * 1e3,
+#   unit = "mt",
+#   uncertainty = fimsdat$cattot_log_sd[1]
+# )
+# index2 <- data.frame(
+#   type = "index",
+#   name = "survey2",
+#   age = NA,
+#   datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"),
+#   dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"),
+#   value = ifelse(ind2 > 0, ind2 * 1e9, ind2),
+#   unit = "",
+#   uncertainty = CV2
+# )
+# index3 <- data.frame(
+#   type = "index",
+#   name = "survey3",
+#   age = NA,
+#   datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"),
+#   dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"),
+#   value = ifelse(ind3 > 0, ind3 * 1e9, ind3),
+#   unit = "",
+#   uncertainty = CV3
+# )
+# index6 <- data.frame(
+#   type = "index",
+#   name = "survey6",
+#   age = NA,
+#   datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"),
+#   dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"),
+#   value = ifelse(ind6 > 0, ind6 * 1e9, ind6),
+#   unit = "",
+#   uncertainty = CV6
+# )
+# ## these have -999 for missing data years
+# catchage <- data.frame(
+#   type = "age",
+#   name = "fleet1",
+#   age = rep(seq(1, nages), nyears),
+#   datestart = rep(paste0(
+#     seq(fimsdat$styr, fimsdat$endyr), "-01-01"
+#   ), each = nages),
+#   dateend = rep(paste0(
+#     seq(fimsdat$styr, fimsdat$endyr), "-12-31"
+#   ), each = nages),
+#   value = as.numeric(t(caa)),
+#   unit = "",
+#   uncertainty = rep(Ncaa, each = nages)
+# )
+# indexage2 <- data.frame(
+#   type = "age",
+#   name = "survey2",
+#   age = rep(seq(1, nages), nyears),
+#   datestart = rep(paste0(
+#     seq(fimsdat$styr, fimsdat$endyr), "-01-01"
+#   ), each = nages),
+#   dateend = rep(paste0(
+#     seq(fimsdat$styr, fimsdat$endyr), "-12-31"
+#   ), each = nages),
+#   value = as.numeric(t(paa2)),
+#   unit = "",
+#   uncertainty = rep(Npaa2, each = nages)
+# )
+# indexage3 <- data.frame(
+#   type = "age",
+#   name = "survey3",
+#   age = rep(seq(1, nages), nyears),
+#   datestart = rep(paste0(
+#     seq(fimsdat$styr, fimsdat$endyr), "-01-01"
+#   ), each = nages),
+#   dateend = rep(paste0(
+#     seq(fimsdat$styr, fimsdat$endyr), "-12-31"
+#   ), each = nages),
+#   value = as.numeric(t(paa3)),
+#   unit = "",
+#   uncertainty = rep(Npaa3, each = nages)
+# )
+# indexage6 <- data.frame(
+#   type = "age",
+#   name = "survey6",
+#   age = rep(seq(1, nages), nyears),
+#   datestart = rep(paste0(
+#     seq(fimsdat$styr, fimsdat$endyr), "-01-01"
+#   ), each = nages),
+#   dateend = rep(paste0(
+#     seq(fimsdat$styr, fimsdat$endyr), "-12-31"
+#   ), each = nages),
+#   value = as.numeric(t(paa6)),
+#   unit = "",
+#   uncertainty = rep(Npaa6, each = nages)
+# )
+# indexage <- rbind(indexage2, indexage3, indexage6)
+# index <- rbind(index2, index3, index6)
+# ## indexage=indexage2
+# ## index=index2
+# res <- rbind(res, landings, index, catchage, indexage)
+# ## rm(landings, index, catchage, indexage)
 
+estimate_fish_selex <- TRUE
+estimate_survey_selex <- TRUE
+estimate_q2 <- TRUE
+estimate_q3 <- TRUE
+estimate_q6 <- TRUE
+estimate_F <- TRUE
+estimate_recdevs <- TRUE
 
-## build a FIMS and PK data set that match
-##  need to fill missing years with -999 so it's ignored in FIMS
-ind2 <- 0 * pkfit0$rep$Eindxsurv2 - 999
-ind2[which(years %in% fimsdat$srvyrs2)] <- fimsdat$indxsurv2
-CV2 <- rep(1, length = nyears) # actually SE in log space
-CV2[which(years %in% fimsdat$srvyrs2)] <- fimsdat$indxsurv_log_sd2
-paa2 <- pkfit0$rep$Esrvp2 * 0 - 999
-paa2[which(years %in% fimsdat$srv_acyrs2), ] <- fimsdat$srvp2
-Npaa2 <- rep(1, nyears)
-Npaa2[which(years %in% fimsdat$srv_acyrs2)] <- fimsdat$multN_srv2
-
-ind3 <- 0 * pkfit0$rep$Eindxsurv3 - 999
-ind3[which(years %in% fimsdat$srvyrs3)] <- fimsdat$indxsurv3
-CV3 <- rep(1, length = nyears) # actually SE in log space
-CV3[which(years %in% fimsdat$srvyrs3)] <- fimsdat$indxsurv_log_sd3
-paa3 <- pkfit0$rep$Esrvp3 * 0 - 999
-paa3[which(years %in% fimsdat$srv_acyrs3), ] <- fimsdat$srvp3
-Npaa3 <- rep(1, nyears)
-Npaa3[which(years %in% fimsdat$srv_acyrs3)] <- fimsdat$multN_srv3
-
-ind6 <- 0 * pkfit0$rep$Eindxsurv6 - 999
-ind6[which(years %in% fimsdat$srvyrs6)] <- fimsdat$indxsurv6
-CV6 <- rep(1, length = nyears) # actually SE in log space
-CV6[which(years %in% fimsdat$srvyrs6)] <- fimsdat$indxsurv_log_sd6
-paa6 <- pkfit0$rep$Esrvp6 * 0 - 999
-paa6[which(years %in% fimsdat$srv_acyrs6), ] <- fimsdat$srvp6
-Npaa6 <- rep(1, nyears)
-Npaa6[which(years %in% fimsdat$srv_acyrs6)] <- fimsdat$multN_srv6
-
-## repeat with fish catch at age, using expected in missing years
-caa <- pkfit0$rep$Ecatp * 0 - 999
-caa[which(years %in% fimsdat$fshyrs), ] <- fimsdat$catp
-Ncaa <- rep(1, nyears)
-Ncaa[which(years %in% fimsdat$fshyrs)] <- fimsdat$multN_fsh
-
-
-
-## put into fims friendly form
-res <- data.frame(
-  type = character(),
-  name = character(),
-  age = integer(),
-  datestart = character(),
-  dateend = character(),
-  value = double(),
-  unit = character(),
-  uncertainty = double()
-)
-landings <- data.frame(
-  type = "landings",
-  name = "fleet1",
-  age = NA,
-  datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"),
-  dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"),
-  value = as.numeric(fimsdat$cattot) * 1e3,
-  unit = "mt",
-  uncertainty = fimsdat$cattot_log_sd[1]
-)
-index2 <- data.frame(
-  type = "index",
-  name = "survey2",
-  age = NA,
-  datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"),
-  dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"),
-  value = ifelse(ind2 > 0, ind2 * 1e9, ind2),
-  unit = "",
-  uncertainty = CV2
-)
-index3 <- data.frame(
-  type = "index",
-  name = "survey3",
-  age = NA,
-  datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"),
-  dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"),
-  value = ifelse(ind3 > 0, ind3 * 1e9, ind3),
-  unit = "",
-  uncertainty = CV3
-)
-index6 <- data.frame(
-  type = "index",
-  name = "survey6",
-  age = NA,
-  datestart = paste0(seq(fimsdat$styr, fimsdat$endyr), "-01-01"),
-  dateend = paste0(seq(fimsdat$styr, fimsdat$endyr), "-12-31"),
-  value = ifelse(ind6 > 0, ind6 * 1e9, ind6),
-  unit = "",
-  uncertainty = CV6
-)
-## these have -999 for missing data years
-catchage <- data.frame(
-  type = "age",
-  name = "fleet1",
-  age = rep(seq(1, nages), nyears),
-  datestart = rep(paste0(
-    seq(fimsdat$styr, fimsdat$endyr), "-01-01"
-  ), each = nages),
-  dateend = rep(paste0(
-    seq(fimsdat$styr, fimsdat$endyr), "-12-31"
-  ), each = nages),
-  value = as.numeric(t(caa)),
-  unit = "",
-  uncertainty = rep(Ncaa, each = nages)
-)
-indexage2 <- data.frame(
-  type = "age",
-  name = "survey2",
-  age = rep(seq(1, nages), nyears),
-  datestart = rep(paste0(
-    seq(fimsdat$styr, fimsdat$endyr), "-01-01"
-  ), each = nages),
-  dateend = rep(paste0(
-    seq(fimsdat$styr, fimsdat$endyr), "-12-31"
-  ), each = nages),
-  value = as.numeric(t(paa2)),
-  unit = "",
-  uncertainty = rep(Npaa2, each = nages)
-)
-indexage3 <- data.frame(
-  type = "age",
-  name = "survey3",
-  age = rep(seq(1, nages), nyears),
-  datestart = rep(paste0(
-    seq(fimsdat$styr, fimsdat$endyr), "-01-01"
-  ), each = nages),
-  dateend = rep(paste0(
-    seq(fimsdat$styr, fimsdat$endyr), "-12-31"
-  ), each = nages),
-  value = as.numeric(t(paa3)),
-  unit = "",
-  uncertainty = rep(Npaa3, each = nages)
-)
-indexage6 <- data.frame(
-  type = "age",
-  name = "survey6",
-  age = rep(seq(1, nages), nyears),
-  datestart = rep(paste0(
-    seq(fimsdat$styr, fimsdat$endyr), "-01-01"
-  ), each = nages),
-  dateend = rep(paste0(
-    seq(fimsdat$styr, fimsdat$endyr), "-12-31"
-  ), each = nages),
-  value = as.numeric(t(paa6)),
-  unit = "",
-  uncertainty = rep(Npaa6, each = nages)
-)
-indexage <- rbind(indexage2, indexage3, indexage6)
-index <- rbind(index2, index3, index6)
-## indexage=indexage2
-## index=index2
-res <- rbind(res, landings, index, catchage, indexage)
-## rm(landings, index, catchage, indexage)
-
-
-age_frame <- FIMS::FIMSFrame(res)
+# age_frame <- FIMS::FIMSFrame(res)
 fishery_catch <- FIMS::m_landings(age_frame)
 fishery_agecomp <- FIMS::m_agecomp(age_frame, "fleet1")
 survey_index2 <- FIMS::m_index(age_frame, "survey2")
@@ -209,118 +214,118 @@ fish_fleet$SetSelectivity(fish_selex$get_id())
 
 
 ## Setup survey 2
-survey_fleet_index <- methods::new(Index, nyears)
-survey_age_comp <- methods::new(AgeComp, nyears, nages)
-survey_fleet_index$index_data <- survey_index2
-survey_age_comp$age_comp_data <-
+survey_fleet_index2 <- methods::new(Index, nyears)
+survey_age_comp2 <- methods::new(AgeComp, nyears, nages)
+survey_fleet_index2$index_data <- survey_index2
+survey_age_comp2$age_comp_data <-
   survey_agecomp2 * indexage2$uncertainty
 ## survey selectivity: ascending logistic
 ## methods::show(DoubleLogisticSelectivity)
-survey_selex <- methods::new(DoubleLogisticSelectivity)
-survey_selex$inflection_point_asc$value <- parfinal$inf1_srv2
-survey_selex$inflection_point_asc$is_random_effect <- FALSE
-survey_selex$inflection_point_asc$estimated <- estimate_survey_selex
-survey_selex$slope_asc$value <- exp(parfinal$log_slp1_srv2)
-survey_selex$slope_asc$is_random_effect <- FALSE
-survey_selex$slope_asc$estimated <- estimate_survey_selex
+survey_selex2 <- methods::new(DoubleLogisticSelectivity)
+survey_selex2$inflection_point_asc$value <- parfinal$inf1_srv2
+survey_selex2$inflection_point_asc$is_random_effect <- FALSE
+survey_selex2$inflection_point_asc$estimated <- estimate_survey_selex
+survey_selex2$slope_asc$value <- exp(parfinal$log_slp1_srv2)
+survey_selex2$slope_asc$is_random_effect <- FALSE
+survey_selex2$slope_asc$estimated <- estimate_survey_selex
 ## not estimated to make it ascending only, fix at input values
-survey_selex$inflection_point_desc$value <- parfinal$inf2_srv2
-survey_selex$inflection_point_desc$is_random_effect <- FALSE
-survey_selex$inflection_point_desc$estimated <- FALSE
-survey_selex$slope_desc$value <- exp(parfinal$log_slp2_srv2)
-survey_selex$slope_desc$is_random_effect <- FALSE
-survey_selex$slope_desc$estimated <- FALSE
-survey_fleet <- methods::new(Fleet)
-survey_fleet$is_survey <- TRUE
-survey_fleet$nages <- nages
-survey_fleet$nyears <- nyears
-survey_fleet$estimate_F <- FALSE
-survey_fleet$random_F <- FALSE
-survey_fleet$log_q <- parfinal$log_q2_mean
-survey_fleet$estimate_q <- estimate_q2
-survey_fleet$random_q <- FALSE
+survey_selex2$inflection_point_desc$value <- parfinal$inf2_srv2
+survey_selex2$inflection_point_desc$is_random_effect <- FALSE
+survey_selex2$inflection_point_desc$estimated <- FALSE
+survey_selex2$slope_desc$value <- exp(parfinal$log_slp2_srv2)
+survey_selex2$slope_desc$is_random_effect <- FALSE
+survey_selex2$slope_desc$estimated <- FALSE
+survey_fleet2 <- methods::new(Fleet)
+survey_fleet2$is_survey <- TRUE
+survey_fleet2$nages <- nages
+survey_fleet2$nyears <- nyears
+survey_fleet2$estimate_F <- FALSE
+survey_fleet2$random_F <- FALSE
+survey_fleet2$log_q <- parfinal$log_q2_mean
+survey_fleet2$estimate_q <- estimate_q2
+survey_fleet2$random_q <- FALSE
 # sd = sqrt(log(cv^2 + 1)), sd is log transformed
-survey_fleet$log_obs_error <- log(index2$uncertainty)
+survey_fleet2$log_obs_error <- log(index2$uncertainty)
 ## survey_fleet$log_obs_error$estimated <- FALSE
-survey_fleet$SetAgeCompLikelihood(1)
-survey_fleet$SetIndexLikelihood(1)
-survey_fleet$SetSelectivity(survey_selex$get_id())
-survey_fleet$SetObservedIndexData(survey_fleet_index$get_id())
-survey_fleet$SetObservedAgeCompData(survey_age_comp$get_id())
+survey_fleet2$SetAgeCompLikelihood(1)
+survey_fleet2$SetIndexLikelihood(1)
+survey_fleet2$SetSelectivity(survey_selex2$get_id())
+survey_fleet2$SetObservedIndexData(survey_fleet_index2$get_id())
+survey_fleet2$SetObservedAgeCompData(survey_age_comp2$get_id())
 
 ## Setup survey 3
-survey_fleet_index <- methods::new(Index, nyears)
-survey_age_comp <- methods::new(AgeComp, nyears, nages)
-survey_fleet_index$index_data <- survey_index3
-survey_age_comp$age_comp_data <-
+survey_fleet_index3 <- methods::new(Index, nyears)
+survey_age_comp3 <- methods::new(AgeComp, nyears, nages)
+survey_fleet_index3$index_data <- survey_index3
+survey_age_comp3$age_comp_data <-
   survey_agecomp3 * indexage3$uncertainty
 ## survey selectivity: ascending logistic
 ## methods::show(LogisticSelectivity)
-survey_selex <- methods::new(LogisticSelectivity)
-survey_selex$inflection_point$value <- parfinal$inf1_srv3
-survey_selex$inflection_point$is_random_effect <- FALSE
-survey_selex$inflection_point$estimated <- estimate_survey_selex
-survey_selex$slope$value <- exp(parfinal$log_slp1_srv3)
-survey_selex$slope$is_random_effect <- FALSE
-survey_selex$slope$estimated <- estimate_survey_selex
-survey_fleet <- methods::new(Fleet)
-survey_fleet$is_survey <- TRUE
-survey_fleet$nages <- nages
-survey_fleet$nyears <- nyears
-survey_fleet$estimate_F <- FALSE
-survey_fleet$random_F <- FALSE
-survey_fleet$log_q <- parfinal$log_q3_mean
-survey_fleet$estimate_q <- estimate_q3
-survey_fleet$random_q <- FALSE
+survey_selex3 <- methods::new(LogisticSelectivity)
+survey_selex3$inflection_point$value <- parfinal$inf1_srv3
+survey_selex3$inflection_point$is_random_effect <- FALSE
+survey_selex3$inflection_point$estimated <- estimate_survey_selex
+survey_selex3$slope$value <- exp(parfinal$log_slp1_srv3)
+survey_selex3$slope$is_random_effect <- FALSE
+survey_selex3$slope$estimated <- estimate_survey_selex
+survey_fleet3 <- methods::new(Fleet)
+survey_fleet3$is_survey <- TRUE
+survey_fleet3$nages <- nages
+survey_fleet3$nyears <- nyears
+survey_fleet3$estimate_F <- FALSE
+survey_fleet3$random_F <- FALSE
+survey_fleet3$log_q <- parfinal$log_q3_mean
+survey_fleet3$estimate_q <- estimate_q3
+survey_fleet3$random_q <- FALSE
 # sd = sqrt(log(cv^2 + 1)), sd is log transformed
-survey_fleet$log_obs_error <- log(index3$uncertainty)
+survey_fleet3$log_obs_error <- log(index3$uncertainty)
 ## survey_fleet$log_obs_error$estimated <- FALSE
-survey_fleet$SetAgeCompLikelihood(2)
-survey_fleet$SetIndexLikelihood(2)
-survey_fleet$SetSelectivity(survey_selex$get_id())
-survey_fleet$SetObservedIndexData(survey_fleet_index$get_id())
-survey_fleet$SetObservedAgeCompData(survey_age_comp$get_id())
+survey_fleet3$SetAgeCompLikelihood(2)
+survey_fleet3$SetIndexLikelihood(2)
+survey_fleet3$SetSelectivity(survey_selex3$get_id())
+survey_fleet3$SetObservedIndexData(survey_fleet_index3$get_id())
+survey_fleet3$SetObservedAgeCompData(survey_age_comp3$get_id())
 
 ## Setup survey 6
-survey_fleet_index <- methods::new(Index, nyears)
-survey_age_comp <- methods::new(AgeComp, nyears, nages)
-survey_fleet_index$index_data <- survey_index6
-survey_age_comp$age_comp_data <-
+survey_fleet_index6 <- methods::new(Index, nyears)
+survey_age_comp6 <- methods::new(AgeComp, nyears, nages)
+survey_fleet_index6$index_data <- survey_index6
+survey_age_comp6$age_comp_data <-
   survey_agecomp6 * indexage6$uncertainty
 ## survey selectivity: ascending logistic
 ## methods::show(DoubleLogisticSelectivity)
-survey_selex <- methods::new(DoubleLogisticSelectivity)
-survey_selex$inflection_point_asc$value <- parfinal$inf1_srv6
-survey_selex$inflection_point_asc$is_random_effect <- FALSE
-survey_selex$inflection_point_asc$estimated <- FALSE
-survey_selex$slope_asc$value <- exp(parfinal$log_slp1_srv6)
-survey_selex$slope_asc$is_random_effect <- FALSE
-survey_selex$slope_asc$estimated <- FALSE
+survey_selex6 <- methods::new(DoubleLogisticSelectivity)
+survey_selex6$inflection_point_asc$value <- parfinal$inf1_srv6
+survey_selex6$inflection_point_asc$is_random_effect <- FALSE
+survey_selex6$inflection_point_asc$estimated <- FALSE
+survey_selex6$slope_asc$value <- exp(parfinal$log_slp1_srv6)
+survey_selex6$slope_asc$is_random_effect <- FALSE
+survey_selex6$slope_asc$estimated <- FALSE
 ## not estimated to make it ascending only, fix at input values
-survey_selex$inflection_point_desc$value <- parfinal$inf2_srv6
-survey_selex$inflection_point_desc$is_random_effect <- FALSE
-survey_selex$inflection_point_desc$estimated <-
+survey_selex6$inflection_point_desc$value <- parfinal$inf2_srv6
+survey_selex6$inflection_point_desc$is_random_effect <- FALSE
+survey_selex6$inflection_point_desc$estimated <-
   estimate_survey_selex
-survey_selex$slope_desc$value <- exp(parfinal$log_slp2_srv6)
-survey_selex$slope_desc$is_random_effect <- FALSE
-survey_selex$slope_desc$estimated <- estimate_survey_selex
-survey_fleet <- methods::new(Fleet)
-survey_fleet$is_survey <- TRUE
-survey_fleet$nages <- nages
-survey_fleet$nyears <- nyears
-survey_fleet$estimate_F <- FALSE
-survey_fleet$random_F <- FALSE
-survey_fleet$log_q <- parfinal$log_q6
-survey_fleet$estimate_q <- estimate_q6
-survey_fleet$random_q <- FALSE
+survey_selex6$slope_desc$value <- exp(parfinal$log_slp2_srv6)
+survey_selex6$slope_desc$is_random_effect <- FALSE
+survey_selex6$slope_desc$estimated <- estimate_survey_selex
+survey_fleet6 <- methods::new(Fleet)
+survey_fleet6$is_survey <- TRUE
+survey_fleet6$nages <- nages
+survey_fleet6$nyears <- nyears
+survey_fleet6$estimate_F <- FALSE
+survey_fleet6$random_F <- FALSE
+survey_fleet6$log_q <- parfinal$log_q6
+survey_fleet6$estimate_q <- estimate_q6
+survey_fleet6$random_q <- FALSE
 # sd = sqrt(log(cv^2 + 1)), sd is log transformed
-survey_fleet$log_obs_error <- log(index6$uncertainty)
+survey_fleet6$log_obs_error <- log(index6$uncertainty)
 ## survey_fleet$log_obs_error$estimated <- FALSE
-survey_fleet$SetAgeCompLikelihood(3)
-survey_fleet$SetIndexLikelihood(3)
-survey_fleet$SetSelectivity(survey_selex$get_id())
-survey_fleet$SetObservedIndexData(survey_fleet_index$get_id())
-survey_fleet$SetObservedAgeCompData(survey_age_comp$get_id())
+survey_fleet6$SetAgeCompLikelihood(3)
+survey_fleet6$SetIndexLikelihood(3)
+survey_fleet6$SetSelectivity(survey_selex6$get_id())
+survey_fleet6$SetObservedIndexData(survey_fleet_index6$get_id())
+survey_fleet6$SetObservedAgeCompData(survey_age_comp6$get_id())
 
 
 
@@ -373,7 +378,7 @@ population$estimate_init_naa <-
   FALSE # TRUE , NOTE: fixing at ASAP estimates to test SSB calculations
 population$nages <- nages
 population$ages <- ages
-population$nfleets <- 2 # 1 fleet and 1 survey
+population$nfleets <- 4 # 1 fleet and 1 survey
 population$nseasons <- nseasons
 population$nyears <- nyears
 ## population$prop_female <- 1.0 # ASAP assumption
diff --git a/content/R/pk_prepare_FIMS_inputs_by_fleet.R b/content/R/pk_prepare_FIMS_inputs_by_fleet.R
index d6ef473..0591b64 100644
--- a/content/R/pk_prepare_FIMS_inputs_by_fleet.R
+++ b/content/R/pk_prepare_FIMS_inputs_by_fleet.R
@@ -1,218 +1,219 @@
-set_selectivity <- function(form, initial_values, is_estimated, is_random_effect) {
-  # Check if the current form is "LogisticSelectivity"
-  if (form == "LogisticSelectivity") {
-    # Create a new LogisticSelectivity object
-    selectivity <- methods::new(LogisticSelectivity)
+# Helper functions -----------------------------------------------
 
-    # Set initial parameter values
-    selectivity$inflection_point$value <- initial_values[1]
-    selectivity$slope$value <- initial_values[2]
-
-    # Set whether each parameter is estimated
-    selectivity$inflection_point$estimated <- is_estimated[1]
-    selectivity$slope$estimated <- is_estimated[2]
+# This function sets the initial parameter values, estimated flags, and random effect flags
+# for a logistic form (e.g., selectivity module or maturity module) within a module.
+use_logistic_form <- function(module, initial_values,
+                              is_estimated, is_random_effect) {
+  # Set initial parameter values for each parameter
+  module$inflection_point$value <- initial_values$inflection_point
+  module$slope$value <- initial_values$slope
 
-    # Set whether each parameter has random effects
-    selectivity$inflection_point$is_random_effect <- is_random_effect[1]
-    selectivity$slope$is_random_effect <- is_random_effect[2]
-  }
+  # Set whether each parameter is estimated
+  module$inflection_point$estimated <- is_estimated$inflection_point
+  module$slope$estimated <- is_estimated$slope
 
-  # Check if the current form is "DoubleLogisticSelectivity"
-  if (form == "DoubleLogisticSelectivity") {
-    # Create a new DoubleLogisticSelectivity object
-    selectivity <- methods::new(DoubleLogisticSelectivity)
+  # Set whether each parameter has random effects
+  module$inflection_point$is_random_effect <- is_random_effect$inflection_point
+  module$slope$is_random_effect <- is_random_effect$slope
 
-    # Set initial parameter values for both ascending and descending parts
-    selectivity$inflection_point_asc$value <- initial_values[1]
-    selectivity$slope_asc$value <- initial_values[2]
-    selectivity$inflection_point_desc$value <- initial_values[3]
-    selectivity$slope_desc$value <- initial_values[4]
+  return(module)
+}
 
-    # Set whether each parameter is estimated
-    selectivity$inflection_point_asc$estimated <- is_estimated[1]
-    selectivity$slope_asc$estimated <- is_estimated[2]
-    selectivity$inflection_point_desc$estimated <- is_estimated[3]
-    selectivity$slope_desc$estimated <- is_estimated[4]
+# This function sets the initial parameter values, estimated flags, and random effect flags
+# for a double logistic form within a module.
+use_double_logistic_form <- function(module, initial_values,
+                                     is_estimated, is_random_effect) {
+  # Set initial parameter values for both ascending and descending parts
+  module$inflection_point_asc$value <- initial_values$inflection_point_asc
+  module$slope_asc$value <- initial_values$slope_asc
+  module$inflection_point_desc$value <- initial_values$inflection_point_desc
+  module$slope_desc$value <- initial_values$slope_desc
 
-    # Set whether each parameter has random effects
-    selectivity$inflection_point_asc$is_random_effect <- is_random_effect[1]
-    selectivity$slope_asc$is_random_effect <- is_random_effect[2]
-    selectivity$inflection_point_desc$is_random_effect <- is_random_effect[3]
-    selectivity$slope_desc$is_random_effect <- is_random_effect[4]
-  }
-  return(selectivity)
+  # Set whether each parameter is estimated
+  module$inflection_point_asc$estimated <- is_estimated$inflection_point_asc
+  module$slope_asc$estimated <- is_estimated$slope_asc
+  module$inflection_point_desc$estimated <- is_estimated$inflection_point_desc
+  module$slope_desc$estimated <- is_estimated$slope_desc
+
+  # Set whether each parameter has random effects
+  module$inflection_point_asc$is_random_effect <- is_random_effect$inflection_point_asc
+  module$slope_asc$is_random_effect <- is_random_effect$slope_asc
+  module$inflection_point_desc$is_random_effect <- is_random_effect$inflection_point_desc
+  module$slope_desc$is_random_effect <- is_random_effect$slope_desc
+
+  return(module)
 }
 
-set_fleet <- function(data = age_frame, is_estimated_obs_error, selectivity_ctl,
-                      selectivity_module_list, Fmort_ctl, catchability_ctl) {
-  fishing_fleet_names <- dplyr::filter(
-    .data = as.data.frame(data@data),
-    type == "landings"
-  ) |>
-    dplyr::distinct(name) |>
-    dplyr::pull(name)
-  survey_fleet_names <- dplyr::filter(
-    .data = as.data.frame(data@data),
-    type == "index"
-  ) |>
-    dplyr::distinct(name) |>
-    dplyr::pull(name)
-
-  fleet_names <- c(fishing_fleet_names, survey_fleet_names)
-  index_types <- c(
-    rep("landings", length(fishing_fleet_names)),
-    rep("index", length(survey_fleet_names))
-  )
-  nfleets <- length(fleet_names)
+# This function sets the fleet data, including selectivity, fishing mortality,
+# catchability, and observational error for each fleet.
+set_fleet <- function(fleet_data, population_data, is_estimated_obs_error, selectivity,
+                      Fmort, catchability) {
+  nfleets <- population_data$nfleets
 
-  nyears <- data@n_years
-  nages <- data@n_ages
+  fishing_fleet_names <- population_data$fishing_fleet_names
+  survey_fleet_names <- population_data$survey_fleet_names
+  fleet_names <- population_data$fleet_names
 
   # Initialize an empty list to store fleet objects
-  index <- age_comp <- age_comp_uncertainty <-
-    selectivity <- fleet <- vector(mode = "list", length = nfleets)
+  index_m <- age_comp_m <- age_comp_uncertainty <-
+    selectivity_m <- fleet_m <- vector(mode = "list", length = nfleets)
 
   for (i in 1:nfleets) {
-    index[[i]] <- methods::new(Index, nyears)
+    index_m[[i]] <- methods::new(Index, nyears)
 
-    if (i %in% seq_along(fishing_fleet_names)) {
-      index[[i]]$index_data <- FIMS::m_landings(data)
+    if (fleet_names[i] %in% fishing_fleet_names) {
+      index_m[[i]]$index_data <- FIMS::m_landings(fleet_data)
     } else {
-      index[[i]]$index_data <- FIMS::m_index(data, fleet_names[i])
+      index_m[[i]]$index_data <- FIMS::m_index(fleet_data, fleet_names[i])
     }
 
     age_comp_uncertainty[[i]] <- dplyr::filter(
-      .data = as.data.frame(data@data),
+      .data = as.data.frame(fleet_data@data),
       name == fleet_names[i] & type == "age"
     ) |>
       dplyr::pull(uncertainty)
 
-    age_comp[[i]] <- methods::new(AgeComp, nyears, nages)
-    age_comp$age_comp_data <- FIMS::m_agecomp(data, fleet_names[i]) *
+    age_comp_m[[i]] <- methods::new(AgeComp, nyears, nages)
+    age_comp_m[[i]]$age_comp_data <- FIMS::m_agecomp(fleet_data, fleet_names[i]) *
       age_comp_uncertainty[[i]]
 
-    if (!is.null(selectivity_ctl)) {
-      selectivity[[i]] <- set_selectivity(
-        form = selectivity_ctl$form[[i]],
-        initial_values = selectivity_ctl$initial_values[[i]],
-        is_estimated = selectivity_ctl$is_estimated[[i]],
-        is_random_effect = selectivity_ctl$is_random_effect[[i]]
+    if (selectivity[[i]]$form == "LogisticSelectivity") {
+      selectivity_temp <- methods::new(LogisticSelectivity)
+      selectivity_m[[i]] <- use_logistic_form(
+        module = selectivity_temp,
+        initial_values = selectivity[[i]]$initial_values,
+        is_estimated = selectivity[[i]]$is_estimated,
+        is_random_effect = selectivity[[i]]$is_random_effect
       )
     }
 
-    if (!is.null(selectivity_module_list)) {
-      selectivity[[i]] <- selectivity_module_list[[i]]
+    if (selectivity[[i]]$form == "DoubleLogisticSelectivity") {
+      selectivity_temp <- methods::new(DoubleLogisticSelectivity)
+      selectivity_m[[i]] <- use_double_logistic_form(
+        module = selectivity_temp,
+        initial_values = selectivity[[i]]$initial_values,
+        is_estimated = selectivity[[i]]$is_estimated,
+        is_random_effect = selectivity[[i]]$is_random_effect
+      )
     }
 
-    fleet[[i]] <- methods::new(Fleet)
+    fleet_m[[i]] <- methods::new(Fleet)
 
     # Set nyears and nages
-    fleet[[i]]$nages <- nages
-    fleet[[i]]$nyears <- nyears
+    fleet_m[[i]]$nages <- population_data$nages
+    fleet_m[[i]]$nyears <- population_data$nyears
 
-    fleet[[i]]$log_obs_error <- dplyr::filter(
-      .data = as.data.frame(data@data),
-      name == fleet_names[i] & type == index_types[i]
+    obs_error <- dplyr::filter(
+      .data = as.data.frame(fleet_data@data),
+      name == fleet_names[i] & type == population_data$index_types[i]
     ) |>
       dplyr::pull(uncertainty)
+    fleet_m[[i]]$log_obs_error <- log(obs_error)
 
-    fleet[[i]]$estimate_obs_error <- is_estimated_obs_error[[i]]
-
-    fleet[[i]]$is_survey <- Fmort_ctl$is_survey[[i]]
-    fleet[[i]]$estimate_F <- Fmort_ctl$estimate_F[[i]]
-    fleet[[i]]$random_F <- Fmort_ctl$random_F[[i]]
+    fleet_m[[i]]$estimate_obs_error <- is_estimated_obs_error[[i]]
 
+    fleet_m[[i]]$is_survey <- Fmort_ctl$is_survey[[i]]
+    fleet_m[[i]]$log_Fmort <- Fmort_ctl$initial_values[[i]]
+    fleet_m[[i]]$estimate_F <- Fmort_ctl$estimate_F[[i]]
+    fleet_m[[i]]$random_F <- Fmort_ctl$random_F[[i]]
 
-    fleet[[i]]$log_q <- catchability_ctl$log_q[[i]]
-    fleet[[i]]$estimate_q <- catchability_ctl$estimate_q[[i]]
-    fleet[[i]]$random_q <- catchability_ctl$random_q[[i]]
+    fleet_m[[i]]$log_q <- catchability_ctl$log_q[[i]]
+    fleet_m[[i]]$estimate_q <- catchability_ctl$estimate_q[[i]]
+    fleet_m[[i]]$random_q <- catchability_ctl$random_q[[i]]
 
-    fleet[[i]]$SetAgeCompLikelihood(i)
-    fleet[[i]]$SetIndexLikelihood(i)
-    fleet[[i]]$SetSelectivity(selectivity[[i]]$get_id())
-    fleet[[i]]$SetObservedIndexData(index[[i]]$get_id())
-    fleet[[i]]$SetObservedAgeCompData(age_comp[[i]]$get_id())
+    fleet_m[[i]]$SetAgeCompLikelihood(i)
+    fleet_m[[i]]$SetIndexLikelihood(i)
+    fleet_m[[i]]$SetSelectivity(selectivity_m[[i]]$get_id())
+    fleet_m[[i]]$SetObservedIndexData(index_m[[i]]$get_id())
+    fleet_m[[i]]$SetObservedAgeCompData(age_comp_m[[i]]$get_id())
   }
-  return(fleet)
+  return(fleet_m)
 }
 
+# This function sets the recruitment parameters, including initial values,
+# estimated flags, and random effect flags, for a specified recruitment form.
 set_recruitment <- function(form, initial_values, is_estimated, is_random_effect) {
   if (form == "BevertonHoltRecruitment") {
-    recruitment <- methods::new(BevertonHoltRecruitment)
+    recruitment_m <- methods::new(BevertonHoltRecruitment)
 
     # Set initial parameter values
-    recruitment$log_rzero$value <- initial_values$log_rzero
-    recruitment$logit_steep$value <- initial_values$logit_steep
-    recruitment$log_sigma_recruit$value <- initial_values$log_sigma_recruit
-    recruitment$log_devs <- initial_values$log_devs
+    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$log_rzero$estimated <- is_estimated$log_rzero
-    recruitment$logit_steep$estimated <- is_estimated$logit_steep
-    recruitment$log_sigma_recruit$estimated <- is_estimated$log_sigma_recruit
-    recruitment$estimate_log_devs <- is_estimated$log_devs
+    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$log_rzero$is_random_effect <- is_random_effect$log_rzero
-    recruitment$logit_steep$is_random_effect <- is_random_effect$logit_steep
-    recruitment$log_sigma_recruit$estimated <- is_estimated$log_sigma_recruit
+    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
   }
-  return(recruitment)
+  return(recruitment_m)
 }
 
-set_growth <- function(data, form, initial_values,
+# This function sets the growth parameters, including initial values,
+# estimated flags, and random effect flags, for a specified growth form.
+set_growth <- function(population_data, form, initial_values,
                        is_estimated = NULL, is_random_effect = NULL) {
   if (form == "EWAAgrowth") {
-    growth <- methods::new(EWAAgrowth)
-    growth$ages <- data$ages
-    # NOTE: FIMS currently cannot use matrix of WAA, so have to ensure constant WAA over time in ASAP file for now
-    growth$weights <- initial_values$weights
+    growth_m <- methods::new(EWAAgrowth)
+    growth_m$ages <- population_data$ages
+    growth_m$weights <- initial_values$weights
   }
-  return(growth)
+  return(growth_m)
 }
 
+# This function sets the maturity parameters, including initial values,
+# estimated flags, and random effect flags, for a specified maturity form.
 set_maturity <- function(form, initial_values, is_estimated, is_random_effect) {
   if (form == "LogisticMaturity") {
-    maturity <- methods::new(LogisticMaturity)
-
-    # Set initial parameter values
-    maturity$inflection_point$value <- initial_values$inflection_point
-    maturity$slope$value <- initial_values$slope
-
-    # Set whether each parameter is estimated
-    maturity$inflection_point$estimated <- is_estimated$inflection_point
-    maturity$slope$estimated <- is_estimated$slope
-
-    # Set whether each parameter has random effects
-    maturity$inflection_point$is_random_effect <- is_random_effect$inflection_point
-    maturity$slope$is_random_effect <- is_random_effect$slope
+    maturity_temp <- methods::new(LogisticMaturity)
+    maturity_m <- use_logistic_form(
+      module = maturity_temp,
+      initial_values = initial_values,
+      is_estimated = is_estimated,
+      is_random_effect = is_random_effect
+    )
   }
 
-  return(maturity)
+  return(maturity_m)
 }
 
-set_population <- function(data, initial_values, is_estimated,
+# This function sets up the population parameters, including initial values,
+# estimated flags, maturity, growth, and recruitment input. Users can use
+# the set_maturity(), set_growth(), and set_recruitment() functions to set up modules
+# before passing them into set_population(). See the example of using set_growth()
+# for setting a module and passing it into set_population().
+set_population <- function(population_data, initial_values, is_estimated,
                            maturity, growth, recruitment) {
-  population <- methods::new(Population)
+  # Create a new population object
+  population_m <- methods::new(Population)
 
-  population$nages <- data$nages
-  population$ages <- data$ages
-  population$nfleets <- data$nfleets
-  population$nseasons <- data$nseasons
-  population$nyears <- data$nyears
+  # Set basic population data
+  population_m$nages <- population_data$nages
+  population_m$ages <- population_data$ages
+  population_m$nfleets <- population_data$nfleets
+  population_m$nseasons <- population_data$nseasons
+  population_m$nyears <- population_data$nyears
 
-  # Set initial parameter values
-  population$log_M <- initial_values$log_M
-  population$log_init_naa <- initial_values$log_init_naa
+  # Set initial values for log_M and log_init_naa
+  population_m$log_M <- initial_values$log_M
+  population_m$log_init_naa <- initial_values$log_init_naa
 
-  # Set whether each parameter is estimated
-  population$estimate_M <- is_estimated$log_M
-  population$estimate_init_naa <- is_estimated$log_init_naa
+  # Set flags indicating which parameters should be estimated
+  population_m$estimate_M <- is_estimated$log_M
+  population_m$estimate_init_naa <- is_estimated$log_init_naa
 
+  # Set up the maturity module
   if (typeof(maturity) == "S4") {
-    maturity_module <- maturity
+    maturity_m <- maturity
   } else {
-    maturity_module <- set_maturity(
+    maturity_m <- set_maturity(
       form = maturity$form,
       initial_values = maturity$initial_values,
       is_estimated = maturity$is_estimated,
@@ -220,20 +221,22 @@ set_population <- function(data, initial_values, is_estimated,
     )
   }
 
+  # Set up the growth module
   if (typeof(growth) == "S4") {
-    growth_module <- growth
+    growth_m <- growth
   } else {
-    growth_module <- set_growth(
-      data = data,
+    growth_m <- set_growth(
+      population_data = population_data,
       form = growth$form,
       initial_values = growth$initial_values
     )
   }
 
+  # Set up the recruitment module
   if (typeof(recruitment) == "S4") {
-    recruitment_module <- recruitment
+    recruitment_m <- recruitment
   } else {
-    recruitment_module <- set_recruitment(
+    recruitment_m <- set_recruitment(
       form = recruitment$form,
       initial_values = recruitment$initial_values,
       is_estimated = recruitment$is_estimated,
@@ -241,9 +244,229 @@ set_population <- function(data, initial_values, is_estimated,
     )
   }
 
-  population$SetMaturity(maturity_module$get_id())
-  population$SetGrowth(growth_module$get_id())
-  population$SetRecruitment(recruitment_module$get_id())
+  # Assign the maturity, growth, and recruitment modules to the population
+  population_m$SetMaturity(maturity_m$get_id())
+  population_m$SetGrowth(growth_m$get_id())
+  population_m$SetRecruitment(recruitment_m$get_id())
 
-  return(population)
+  return(population_m)
 }
+
+# Examples for using helper functions -------------------------------------
+
+# set up fleet data
+fishing_fleet_names <- dplyr::filter(
+  .data = as.data.frame(age_frame@data),
+  type == "landings"
+) |>
+  dplyr::distinct(name) |>
+  dplyr::pull(name)
+survey_fleet_names <- dplyr::filter(
+  .data = as.data.frame(age_frame@data),
+  type == "index"
+) |>
+  dplyr::distinct(name) |>
+  dplyr::pull(name)
+
+fleet_names <- c(fishing_fleet_names, survey_fleet_names)
+index_types <- c(
+  rep("landings", length(fishing_fleet_names)),
+  rep("index", length(survey_fleet_names))
+)
+nfleets <- length(fleet_names)
+
+# set up population data
+population_data <- list(
+  nages = age_frame@n_ages,
+  ages = age_frame@ages,
+  nfleets = nfleets,
+  fishing_fleet_names = fishing_fleet_names,
+  survey_fleet_names = survey_fleet_names,
+  fleet_names = fleet_names,
+  index_types = index_types,
+  nseasons = 1,
+  nyears = age_frame@n_years
+)
+
+# Define whether observation error should be estimated for each fleet
+is_estimated_obs_error <- list(FALSE, FALSE, FALSE, FALSE)
+
+# Define selectivity for fishing fleet and surveys using Double Logistic and Logistic forms
+selectivity_fsh_mean <-
+  selectivity_srv2 <- selectivity_srv6 <-
+  list(
+    form = list("DoubleLogisticSelectivity"),
+    initial_values = list(
+      inflection_point_asc = parfinal$inf1_fsh_mean,
+      slope_asc = exp(parfinal$log_slp1_fsh_mean),
+      inflection_point_desc = parfinal$inf2_fsh_mean,
+      slope_desc = exp(parfinal$log_slp2_fsh_mean)
+    ),
+    is_estimated = list(
+      inflection_point_asc = TRUE,
+      slope_asc = TRUE,
+      inflection_point_desc = TRUE,
+      slope_desc = TRUE
+    ),
+    is_random_effect = list(
+      inflection_point_asc = FALSE,
+      slope_asc = FALSE,
+      inflection_point_desc = FALSE,
+      slope_desc = FALSE
+    )
+  )
+
+# Adjust initial values for selectivity for Survey 2
+selectivity_srv2$initial_values <- list(
+  inflection_point_asc = parfinal$inf1_srv2,
+  slope_asc = exp(parfinal$log_slp1_srv2),
+  inflection_point_desc = parfinal$inf2_srv2,
+  slope_desc = exp(parfinal$log_slp2_srv2)
+)
+
+# Set estimated flags for descending inflection point and slope to FALSE for Survey 2
+selectivity_srv2$is_estimated$inflection_point_desc <-
+  selectivity_srv2$is_estimated$slope_desc <- FALSE
+
+# Adjust initial values for selectivity for Survey 6
+selectivity_srv6$initial_values <- list(
+  inflection_point_asc = parfinal$inf1_srv6,
+  slope_asc = exp(parfinal$log_slp1_srv6),
+  inflection_point_desc = parfinal$inf2_srv6,
+  slope_desc = exp(parfinal$log_slp2_srv6)
+)
+
+# Set estimated flags for ascending inflection point and slope to FALSE for Survey 6
+selectivity_srv6$is_estimated$inflection_point_asc <-
+  selectivity_srv6$is_estimated$slope_asc <- FALSE
+
+# Define selectivity for Survey 3 using Logistic form
+selectivity_srv3 <- list(
+  form = list("LogisticSelectivity"),
+  initial_values = list(
+    inflection_point = parfinal$inf1_srv3,
+    slope = exp(parfinal$log_slp1_srv3)
+  ),
+  is_estimated = list(
+    inflection_point = TRUE,
+    slope = TRUE
+  ),
+  is_random_effect = list(
+    inflection_point = FALSE,
+    slope = FALSE
+  )
+)
+
+# Combine all selectivity controls into a list (order matters!)
+selectivity_ctl <- list(
+  selectivity_fsh_mean,
+  selectivity_srv2,
+  selectivity_srv3,
+  selectivity_srv6
+)
+
+# Define fishing mortality controls
+Fmort_ctl <- list(
+  is_survey = list(FALSE, TRUE, TRUE, TRUE),
+  initial_values = list(log(pkfitfinal$rep$F), 0, 0, 0),
+  estimate_F = list(TRUE, FALSE, FALSE, FALSE),
+  random_F = list(FALSE, FALSE, FALSE, FALSE)
+)
+
+# Define catchability controls
+catchability_ctl <- list(
+  log_q = list(
+    0, parfinal$log_q2_mean,
+    parfinal$log_q3_mean, parfinal$log_q6
+  ),
+  estimate_q = list(FALSE, TRUE, TRUE, TRUE),
+  random_q = list(FALSE, FALSE, FALSE, FALSE)
+)
+
+# Set up the fleet module using the defined controls
+fleet <- set_fleet(
+  fleet_data = age_frame,
+  population_data = population_data,
+  is_estimated_obs_error = is_estimated_obs_error,
+  selectivity = selectivity_ctl,
+  Fmort = Fmort_ctl,
+  catchability = catchability_ctl
+)
+
+# Define recruitment controls using Beverton-Holt form
+recruitment_ctl <- list(
+  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
+  )
+)
+
+# Define growth controls using EWAAgrowth form
+# growth_ctl <- list(
+#   data = population_data,
+#   form = "EWAAgrowth",
+#   initial_values = list(weights = pkinput$dat$wt_srv1[1, ])
+# )
+
+growth_m <- set_growth(
+  population_data = population_data,
+  form = "EWAAgrowth",
+  initial_values = list(weights = pkinput$dat$wt_srv1[1, ])
+)
+
+# Define maturity controls using logistic form
+maturity_ctl <- list(
+  form = "LogisticMaturity",
+  initial_values = list(
+    inflection_point = 4.5,
+    slope = 1.5
+  ),
+  is_estimated = list(
+    inflection_point = FALSE,
+    slope = FALSE
+  ),
+  is_random_effect = list(
+    inflection_point = FALSE,
+    slope = FALSE
+  )
+)
+
+# Define initial values for log_M and log__init_naa
+population_initial_values <- list(
+  log_M = log(as.numeric(t(matrix(
+    rep(pkfitfinal$rep$M, each = population_data$nyears),
+    nrow = population_data$nyears
+  )))),
+  log_init_naa =
+    c(log(pkfitfinal$rep$recruit[1]), log(pkfitfinal$rep$initN)) + log(1e9)
+)
+
+# Define estimated flags for log_M and log__init_naa
+population_is_estimated <- list(
+  log_M = FALSE,
+  log_init_naa = FALSE
+)
+
+# Set up the population module using the defined controls and initial values
+population <- set_population(
+  population_data = population_data,
+  initial_values = population_initial_values,
+  is_estimated = population_is_estimated,
+  maturity = maturity_ctl,
+  growth = growth_m,
+  recruitment = recruitment_ctl
+)
diff --git a/content/R/pk_prepare_FIMS_inputs_by_process.R b/content/R/pk_prepare_FIMS_inputs_by_process.R
index 5509b17..3ad4332 100644
--- a/content/R/pk_prepare_FIMS_inputs_by_process.R
+++ b/content/R/pk_prepare_FIMS_inputs_by_process.R
@@ -4,14 +4,14 @@
 
 
 setup_fleets <- function(fleets,population, catchability, selectivity, data){
-  flts <- list()
+  flts <- fltindex <- fltacomp <- list()
   for(i in 1:length(fleets$names)){
     isfleet <- fleets$type[i] == 'landings'
     flts[[i]] <- methods::new(Fleet)
     flts[[i]]$nages <- length(population$ages)
     flts[[i]]$nyears <- population$nyears
-    fltindex <- methods::new(Index, nyears)
-    fltacomp <- methods::new(AgeComp, flts[[i]]$nyears, flts[[i]]$nages)
+    fltindex[[i]] <- methods::new(Index, nyears)
+    fltacomp[[i]] <- methods::new(AgeComp, flts[[i]]$nyears, flts[[i]]$nages)
     fltdat <- filter(data@data, name==fleets$names[i])
     acomp <- FIMS::m_agecomp(age_frame, fleets$names[i])
     if(isfleet){
@@ -19,21 +19,21 @@ setup_fleets <- function(fleets,population, catchability, selectivity, data){
       flts[[i]]$log_Fmort <- log(fleets$init_Fmort[[i]])
       flts[[i]]$estimate_F <- fleets$estimate_F[i]
       flts[[i]]$random_F <- fleets$reFmort[i]
-      fltindex$index_data <- FIMS::m_landings(data)
+      fltindex[[i]]$index_data <- FIMS::m_landings(data)
       flts[[i]]$log_obs_error <- log(filter(fltdat, type=='landings')$uncertainty)
-      fltacomp$age_comp_data <- acomp * filter(fltdat, type=='age')$uncertainty
+      fltacomp[[i]]$age_comp_data <- acomp * filter(fltdat, type=='age')$uncertainty
     } else {
       flts[[i]]$is_survey <- TRUE
       flts[[i]]$log_q <- catchability$init_pars[i]
       flts[[i]]$estimate_q <- is.na(catchability$fix_pars[i])
       flts[[i]]$random_q <- FALSE# catchability$re[i]!='none'
-      fltindex$index_data <- FIMS::m_index(data, fleets$names[i])
+      fltindex[[i]]$index_data <- FIMS::m_index(data, fleets$names[i])
       flts[[i]]$log_obs_error <- log(filter(fltdat, type=='index')$uncertainty)
-      fltacomp$age_comp_data <- acomp * filter(fltdat, type=='age')$uncertainty
+      fltacomp[[i]]$age_comp_data <- acomp * filter(fltdat, type=='age')$uncertainty
     }
     ## Set Index, AgeComp, and Selectivity using the IDs from the modules defined above
-    flts[[i]]$SetObservedIndexData(fltindex$get_id())
-    flts[[i]]$SetObservedAgeCompData(fltacomp$get_id())
+    flts[[i]]$SetObservedIndexData(fltindex[[i]]$get_id())
+    flts[[i]]$SetObservedAgeCompData(fltacomp[[i]]$get_id())
     flts[[i]]$SetSelectivity(selectivity[[i]]$get_id())
   }
   message("Finished setting up fleets")
diff --git a/content/R/pk_prepare_dat.R b/content/R/pk_prepare_dat.R
index 680f0d4..bd140fe 100644
--- a/content/R/pk_prepare_dat.R
+++ b/content/R/pk_prepare_dat.R
@@ -1,5 +1,3 @@
-
-
 ## build a FIMS and PK data set that match
 
 pkfitfinal <- readRDS(file.path(getwd(), "content", "data_files", "pkfitfinal.RDS"))
@@ -161,20 +159,3 @@ res <- rbind(res, landings, index, catchage, indexage)
 ## rm(landings, index, catchage, indexage)
 
 age_frame <- FIMS::FIMSFrame(res)
-fishery_catch <- FIMS::m_landings(age_frame)
-fishery_agecomp <- FIMS::m_agecomp(age_frame, "fleet1")
-survey_index2 <- FIMS::m_index(age_frame, "survey2")
-survey_agecomp2 <- FIMS::m_agecomp(age_frame, "survey2")
-survey_index3 <- FIMS::m_index(age_frame, "survey3")
-survey_agecomp3 <- FIMS::m_agecomp(age_frame, "survey3")
-survey_index6 <- FIMS::m_index(age_frame, "survey6")
-survey_agecomp6 <- FIMS::m_agecomp(age_frame, "survey6")
-# need to think about how to deal with multiple fleets - only using 1 fleeet for now
-fish_index <- methods::new(Index, nyears)
-fish_age_comp <- methods::new(AgeComp, nyears, nages)
-fish_index$index_data <- fishery_catch
-fish_age_comp$age_comp_data <-
-  fishery_agecomp * catchage$uncertainty#rep(Ncaa, each=nages)
-
-
-FIMS::m_index(age_frame, 'survey2')
diff --git a/content/run_pollock_tests.R b/content/run_pollock_tests.R
index e82f952..8eafd66 100644
--- a/content/run_pollock_tests.R
+++ b/content/run_pollock_tests.R
@@ -1,66 +1,88 @@
-
-
+#  Load packages and input data -------------------------------------------
+# Load required packages
 packages <- c("dplyr", "tidyr", "ggplot2", "TMB", "reshape2", "here", "remotes", "lubridate")
 invisible(lapply(packages, library, character.only = TRUE))
-## install.packages("FIMS", repos = c("https://noaa-fims.r-universe.dev", "https://cloud.r-project.org"))
+
+# Install the FIMS package from specific repositories
+# install.packages("FIMS", repos = c("https://noaa-fims.r-universe.dev", "https://cloud.r-project.org"))
 library(FIMS)
+
+# Define the years and ages for the assessment
 years <- 1970:2023
 nyears <- length(years)
 nseasons <- 1
 nages <- 10
 ages <- 1:nages
-## nfleets <- 1
-## This will fit the models bridging to FIMS (simplifying)
-## source("fit_bridge_models.R")
-## compare changes to model
-# source("R/pk_prepare_dat.R")
-source(file.path(getwd(), "content", "R", "pk_prepare_dat.R"))
-## Some global settings which I Think we can ignore for now
-estimate_fish_selex <- TRUE
-estimate_survey_selex <- TRUE
-estimate_q2 <- TRUE
-estimate_q3 <- TRUE
-estimate_q6 <- TRUE
-estimate_F <- TRUE
-estimate_recdevs <- TRUE
-
 
+# Source the script to prepare input data
+source(file.path(getwd(), "content", "R", "pk_prepare_dat.R"))
 
+# Run FIMS without helper functions ---------------------------------------
 clear()
 clear_logs()
-# source("R/pk_prepare_FIMS_inputs.R")
+
+# Source the script to prepare FIMS inputs without using helper functions
 source(file.path(getwd(), "content", "R", "pk_prepare_FIMS_inputs.R"))
-## make FIMS model
+
+# Create the TMB model and generate the report
 success <- CreateTMBModel()
 parameters <- list(p = get_fixed())
-obj1 <- MakeADFun(data = list(), parameters, DLL = "FIMS", silent = TRUE)
+obj1 <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS", silent = TRUE)
 rep1 <- obj1$report()
 
-
-## Organize input lists by process
+# Run FIMS with helper functions (by process) -------------------------------
 clear()
 clear_logs()
-# source("R/pk_prepare_FIMS_inputs_by_process.R")
+
 source(file.path(getwd(), "content", "R", "pk_prepare_FIMS_inputs_by_process.R"))
-## make FIMS model
+
 success <- CreateTMBModel()
 parameters <- list(p = get_fixed())
 obj2 <- MakeADFun(data = list(), parameters, DLL = "FIMS", silent = TRUE)
-## report values for the two models
+
 rep2 <- obj2$report()
 
-## Test that the two models are identical
-all.equal(rep1,rep2)
+# Run FIMS with helper functions (by fleet) ---------------------------------
+clear()
+clear_logs()
+source(file.path(getwd(), "content", "R", "pk_prepare_FIMS_inputs_by_fleet.R"))
+
+success <- CreateTMBModel()
+parameters <- list(p = get_fixed())
+obj3 <- MakeADFun(data = list(), parameters, DLL = "FIMS", silent = TRUE)
+
+rep3 <- obj3$report()
+
+
+# Compare the reports from the three different runs -----------------------
+# Check if reports from the first two runs are identical (original v.s. by process)
+all.equal(rep1, rep2)
+# Check if reports from the first and third runs are identical (original v.s. by fleet)
+all.equal(rep1, rep3)
 
-cbind(sort(obj1$par), sort(obj2$par))
+# Print the sorted parameter values from the three runs
+print(cbind(sort(obj1$par), sort(obj2$par), sort(obj3$par)))
 
 clear()
 clear_logs()
 
-## Parameters are in a different order but seem to match
-opt1 <- with(obj1, nlminb(par,fn,gr))
-opt2 <- with(obj2, nlminb(par,fn,gr))
+# Optimize the models and compare the optimized parameters
+opt1 <- with(obj1, nlminb(par, fn, gr))
+opt2 <- with(obj2, nlminb(par, fn, gr))
+opt3 <- with(obj3, nlminb(par, fn, gr))
+
+# Compare optimized reports
 all.equal(obj1$report(opt1$par), obj2$report(opt2$par))
-sum(obj1$par)-sum(obj2$par)
+all.equal(obj1$report(opt1$par), obj3$report(opt2$par))
+
+# Check the sum of parameters difference
+sum(obj1$par) - sum(obj2$par)
+sum(obj1$par) - sum(obj3$par)
+
+# Print the sorted optimized parameter values from the three runs
+print(cbind(sort(obj1$par), sort(obj2$par), sort(obj3$par)))
+
+clear()
+clear_logs()
 
-cbind(sort(opt1$par)- sort(opt2$par))
+rm(list = ls())
diff --git a/content/run_pollock_tests_by_fleet.R b/content/run_pollock_tests_by_fleet.R
index f51c3b2..e479404 100644
--- a/content/run_pollock_tests_by_fleet.R
+++ b/content/run_pollock_tests_by_fleet.R
@@ -1,184 +1,20 @@
 ## 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"))
-
+# 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()
 
-is_estimated_obs_error <- list(FALSE, FALSE, FALSE, FALSE)
-selectivity_ctl <- list(
-  form = list(
-    "DoubleLogisticSelectivity",
-    "DoubleLogisticSelectivity",
-    "LogisticSelectivity",
-    "DoubleLogisticSelectivity"
-  ),
-  initial_values = list(
-    c(
-      parfinal$inf1_fsh_mean, exp(parfinal$log_slp1_fsh_mean),
-      parfinal$inf2_fsh_mean, exp(parfinal$log_slp2_fsh_mean)
-    ),
-    c(
-      parfinal$inf1_srv2, exp(parfinal$log_slp1_srv2),
-      parfinal$inf2_srv2, exp(parfinal$log_slp2_srv2)
-    ),
-    c(parfinal$inf1_srv3, exp(parfinal$log_slp1_srv3)),
-    c(
-      parfinal$inf1_srv6, exp(parfinal$log_slp1_srv6),
-      parfinal$inf2_srv6, exp(parfinal$log_slp2_srv6)
-    )
-  ),
-  is_estimated = list(
-    rep(TRUE, 4),
-    rep(TRUE, 4),
-    rep(TRUE, 2),
-    rep(TRUE, 4)
-  ),
-  is_random_effect = list(
-    rep(FALSE, 4),
-    rep(FALSE, 4),
-    rep(FALSE, 2),
-    rep(FALSE, 4)
-  )
-)
-
-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)
-)
-
-fleet <- set_fleet(
-  data = age_frame,
-  is_estimated_obs_error = is_estimated_obs_error,
-  selectivity_ctl = selectivity_ctl,
-  selectivity_module_list = NULL,
-  Fmort_ctl = Fmort_ctl,
-  catchability_ctl = catchability_ctl
-)
-
-# Population module
-# population
-fishing_fleet_names <- dplyr::filter(
-  .data = as.data.frame(age_frame@data),
-  type == "landings"
-) |>
-  dplyr::distinct(name) |>
-  dplyr::pull(name)
-survey_fleet_names <- dplyr::filter(
-  .data = as.data.frame(age_frame@data),
-  type == "index"
-) |>
-  dplyr::distinct(name) |>
-  dplyr::pull(name)
-
-fleet_names <- c(fishing_fleet_names, survey_fleet_names)
-index_types <- c(
-  rep("landings", length(fishing_fleet_names)),
-  rep("index", length(survey_fleet_names))
-)
-nfleets <- length(fleet_names)
-
-population_data <- list(
-  nages = age_frame@n_ages,
-  ages = age_frame@ages,
-  nfleets = nfleets,
-  nseasons = 1,
-  nyears = age_frame@n_years
-)
-
-# recruitment
-recruitment_ctl <- list(
-  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_ctl <- list(
-  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_ctl <- list(
-  form = "LogisticMaturity",
-  initial_values = list(
-    inflection_point = 4.5,
-    slope = 1.5
-  ),
-  is_estimated = list(
-    inflection_point = FALSE,
-    slope = FALSE
-  ),
-  is_random_effect = list(
-    inflection_point = FALSE,
-    slope = FALSE
-  )
-)
-
-population_initial_values <- list(
-  log_M = log(as.numeric(t(matrix(
-    rep(pkfitfinal$rep$M, each = nyears),
-    nrow = nyears
-  )))),
-  log_init_naa =
-    c(log(pkfitfinal$rep$recruit[1]), log(pkfitfinal$rep$initN)) + log(1e9)
-)
-
-population_is_estimated <- list(
-  log_M = FALSE,
-  log_init_naa = FALSE
-)
-population <- set_population(
-  data = population_data,
-  initial_values = population_initial_values,
-  is_estimated = population_is_estimated,
-  maturity = maturity_ctl,
-  growth = growth_ctl,
-  recruitment = recruitment_ctl
-)
 
 ## make FIMS model
 success <- CreateTMBModel()
@@ -187,12 +23,38 @@ 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",