From cad92a3ea5daf444c2b3d80f5f07ec066e76a6b1 Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA <85530309+Andrea-Havron-NOAA@users.noreply.github.com> Date: Tue, 31 Oct 2023 19:09:59 +0000 Subject: [PATCH 1/4] add documentation Co-authored-by: KyleShertzer-NOAA Co-authored-by: Ian Taylor --- tests/testthat/test-fims-estimation.R | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/tests/testthat/test-fims-estimation.R b/tests/testthat/test-fims-estimation.R index 20344b303..ca7b6adbe 100644 --- a/tests/testthat/test-fims-estimation.R +++ b/tests/testthat/test-fims-estimation.R @@ -24,20 +24,37 @@ setwd(working_dir) on.exit(setwd(working_dir), add = TRUE) # Set-up Rcpp modules and fix parameters to "true" setup_fims <- function(om_input, om_output, em_input) { + # create new R environment in which to conduct tests test_env <- new.env(parent = emptyenv()) + # initialize FIMS model test_env$fims <- Rcpp::Module("fims", PACKAGE = "FIMS") # Recruitment - test_env$recruitment <- new(test_env$fims$BevertonHoltRecruitment) - # logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is taken - # before the likelihood calculation + # create new module in the recruitment class (specifically Beverton-Holt, + # when there are other option, this would be where the option would be chosen) + test_env$recruitment <- methods::new(test_env$fims$BevertonHoltRecruitment) + + # NOTE: in first set of parameters below (for recruitment), + # $is_random_effect (default is FALSE) and $estimated (default is FALSE) + # are defined even if they match the defaults in order to provide an example + # of how that is done. Other sections of the code below leave defaults in + # place as appropriate. + + # set up logR_sd + # logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is + # taken before the likelihood calculation test_env$recruitment$log_sigma_recruit$value <- log(om_input$logR_sd) + test_env$recruitment$log_sigma$is_random_effect <- FALSE + test_env$recruitment$log_sigma$estimated <- FALSE + # set up log_rzero (equilibrium recruitment) test_env$recruitment$log_rzero$value <- log(om_input$R0) test_env$recruitment$log_rzero$is_random_effect <- FALSE test_env$recruitment$log_rzero$estimated <- TRUE + # set up logit_steep test_env$recruitment$logit_steep$value <- -log(1.0 - om_input$h) + log(om_input$h - 0.2) test_env$recruitment$logit_steep$is_random_effect <- FALSE test_env$recruitment$logit_steep$estimated <- FALSE + # turn on estimation of deviations test_env$recruitment$estimate_deviations <- TRUE # recruit deviations should enter the model in normal space. # The log is taken in the likelihood calculations From 9c22e6803c89dfc2bbd791a7adbc280ab3ca5b21 Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA Date: Tue, 7 Nov 2023 16:05:59 -0800 Subject: [PATCH 2/4] add more documentation Co-authored-by: Ian Taylor Co-authored-by: KyleShertzer-NOAA --- tests/testthat/test-fims-estimation.R | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/tests/testthat/test-fims-estimation.R b/tests/testthat/test-fims-estimation.R index ca7b6adbe..f73813ef2 100644 --- a/tests/testthat/test-fims-estimation.R +++ b/tests/testthat/test-fims-estimation.R @@ -16,6 +16,8 @@ FIMS_C0_estimation <- ASSAMC::save_initial_input( case_name = "FIMS_C0_estimation" ) +# generate om_input, om_output, and em_input +# using function from the model comparison project ASSAMC::run_om(input_list = FIMS_C0_estimation) on.exit(unlink(maindir, recursive = T), add = TRUE) @@ -23,6 +25,7 @@ on.exit(unlink(maindir, recursive = T), add = TRUE) setwd(working_dir) on.exit(setwd(working_dir), add = TRUE) # Set-up Rcpp modules and fix parameters to "true" +# om_input, om_output, and em_input generated from ASSAMC::run_om setup_fims <- function(om_input, om_output, em_input) { # create new R environment in which to conduct tests test_env <- new.env(parent = emptyenv()) @@ -31,7 +34,7 @@ setup_fims <- function(om_input, om_output, em_input) { # Recruitment # create new module in the recruitment class (specifically Beverton-Holt, - # when there are other option, this would be where the option would be chosen) + # when there are other options, this would be where the option would be chosen) test_env$recruitment <- methods::new(test_env$fims$BevertonHoltRecruitment) # NOTE: in first set of parameters below (for recruitment), @@ -46,7 +49,7 @@ setup_fims <- function(om_input, om_output, em_input) { test_env$recruitment$log_sigma_recruit$value <- log(om_input$logR_sd) test_env$recruitment$log_sigma$is_random_effect <- FALSE test_env$recruitment$log_sigma$estimated <- FALSE - # set up log_rzero (equilibrium recruitment) + # set up log_rzero (equilibrium recruitment) test_env$recruitment$log_rzero$value <- log(om_input$R0) test_env$recruitment$log_rzero$is_random_effect <- FALSE test_env$recruitment$log_rzero$estimated <- TRUE @@ -63,11 +66,15 @@ setup_fims <- function(om_input, om_output, em_input) { # Data test_env$catch <- em_input$L.obs$fleet1 + # set fishing fleet catch data, need to set dimensions of data index + # currently FIMS only has a fleet module that takes index for both survey index and fishery catch test_env$fishing_fleet_index <- new(test_env$fims$Index, length(test_env$catch)) test_env$fishing_fleet_index$index_data <- test_env$catch + # set fishing fleet age comp data, need to set dimensions of age comps test_env$fishing_fleet_age_comp <- new(test_env$fims$AgeComp, length(test_env$catch), om_input$nages) test_env$fishing_fleet_age_comp$age_comp_data <- c(t(em_input$L.age.obs$fleet1)) * em_input$n.L$fleet1 + # repeat for surveys test_env$survey_index <- em_input$surveyB.obs$survey1 test_env$survey_fleet_index <- new(test_env$fims$Index, length(test_env$survey_index)) test_env$survey_fleet_index$index_data <- test_env$survey_index @@ -93,8 +100,10 @@ setup_fims <- function(om_input, om_output, em_input) { test_env$fishing_fleet_selectivity <- new(test_env$fims$LogisticSelectivity) test_env$fishing_fleet_selectivity$inflection_point$value <- om_input$sel_fleet$fleet1$A50.sel1 test_env$fishing_fleet_selectivity$inflection_point$is_random_effect <- FALSE + # turn on estimation of inflection_point test_env$fishing_fleet_selectivity$inflection_point$estimated <- TRUE test_env$fishing_fleet_selectivity$slope$value <- om_input$sel_fleet$fleet1$slope.sel1 + # turn on estimation of slope test_env$fishing_fleet_selectivity$slope$is_random_effect <- FALSE test_env$fishing_fleet_selectivity$slope$estimated <- TRUE @@ -109,7 +118,10 @@ setup_fims <- function(om_input, om_output, em_input) { test_env$fishing_fleet$random_q <- FALSE test_env$fishing_fleet$log_obs_error$value <- log(sqrt(log(em_input$cv.L$fleet1^2 + 1))) test_env$fishing_fleet$log_obs_error$estimated <- FALSE - # Need get_id() for setting up observed agecomp and index data? + # Modules are linked together using module IDs + # Each module has a get_id() function that returns the unique ID for that module + # Each fleet uses the module IDs to link up the correct module to the correct fleet + # Note: Likelihoods not yet set up as a stand-alone modules, so no get_id() test_env$fishing_fleet$SetAgeCompLikelihood(1) test_env$fishing_fleet$SetIndexLikelihood(1) test_env$fishing_fleet$SetSelectivity(test_env$fishing_fleet_selectivity$get_id()) @@ -120,9 +132,11 @@ setup_fims <- function(om_input, om_output, em_input) { test_env$survey_fleet_selectivity <- new(test_env$fims$LogisticSelectivity) test_env$survey_fleet_selectivity$inflection_point$value <- om_input$sel_survey$survey1$A50.sel1 test_env$survey_fleet_selectivity$inflection_point$is_random_effect <- FALSE + # turn on estimation of inflection_point test_env$survey_fleet_selectivity$inflection_point$estimated <- TRUE test_env$survey_fleet_selectivity$slope$value <- om_input$sel_survey$survey1$slope.sel1 test_env$survey_fleet_selectivity$slope$is_random_effect <- FALSE + # turn on estimation of slope test_env$survey_fleet_selectivity$slope$estimated <- TRUE test_env$survey_fleet <- new(test_env$fims$Fleet) From 5e9b9455bca0dff992bfeb3ca3281a5e26f01c6f Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA Date: Tue, 7 Nov 2023 16:24:55 -0800 Subject: [PATCH 3/4] fix parameter name --- tests/testthat/test-fims-estimation.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-fims-estimation.R b/tests/testthat/test-fims-estimation.R index f73813ef2..168363671 100644 --- a/tests/testthat/test-fims-estimation.R +++ b/tests/testthat/test-fims-estimation.R @@ -47,8 +47,8 @@ setup_fims <- function(om_input, om_output, em_input) { # logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is # taken before the likelihood calculation test_env$recruitment$log_sigma_recruit$value <- log(om_input$logR_sd) - test_env$recruitment$log_sigma$is_random_effect <- FALSE - test_env$recruitment$log_sigma$estimated <- FALSE + test_env$recruitment$log_sigma_recruit$is_random_effect <- FALSE + test_env$recruitment$log_sigma_recruit$estimated <- FALSE # set up log_rzero (equilibrium recruitment) test_env$recruitment$log_rzero$value <- log(om_input$R0) test_env$recruitment$log_rzero$is_random_effect <- FALSE From 5d431ad4df523334f69aac88ee00614c8435617e Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA <85530309+Andrea-Havron-NOAA@users.noreply.github.com> Date: Thu, 16 Nov 2023 18:01:19 +0000 Subject: [PATCH 4/4] add more documentation and clean up test Co-authored-by: Ian Taylor Co-authored-by: KyleShertzer-NOAA --- tests/testthat/test-fims-estimation.R | 60 +++++++++++++++------------ 1 file changed, 34 insertions(+), 26 deletions(-) diff --git a/tests/testthat/test-fims-estimation.R b/tests/testthat/test-fims-estimation.R index ec76f8e4c..d79a02dce 100644 --- a/tests/testthat/test-fims-estimation.R +++ b/tests/testthat/test-fims-estimation.R @@ -20,7 +20,7 @@ FIMS_C0_estimation <- ASSAMC::save_initial_input( # using function from the model comparison project ASSAMC::run_om(input_list = FIMS_C0_estimation) -on.exit(unlink(maindir, recursive = T), add = TRUE) +on.exit(unlink(maindir, recursive = TRUE), add = TRUE) setwd(working_dir) on.exit(setwd(working_dir), add = TRUE) @@ -144,7 +144,6 @@ setup_fims <- function(om_input, om_output, em_input) { test_env$survey_fleet$is_survey <- TRUE test_env$survey_fleet$nages <- om_input$nages test_env$survey_fleet$nyears <- om_input$nyr - # survey_fleet$log_Fmort <- rep(log(0.0000000000000000000000000001), om_input$nyr) #-Inf? test_env$survey_fleet$estimate_F <- FALSE test_env$survey_fleet$random_F <- FALSE test_env$survey_fleet$log_q <- log(om_output$survey_q$survey1) @@ -160,9 +159,6 @@ setup_fims <- function(om_input, om_output, em_input) { # Population test_env$population <- new(test_env$fims$Population) - # is it a problem these are not Parameters in the Population interface? - # the Parameter class (from rcpp/rcpp_objects/rcpp_interface_base) cannot handle vectors, - # do we need a ParameterVector class? test_env$population$log_M <- rep(log(om_input$M.age[1]), om_input$nyr * om_input$nages) test_env$population$estimate_M <- FALSE test_env$population$log_init_naa <- log(om_output$N.age[1, ]) @@ -172,30 +168,36 @@ setup_fims <- function(om_input, om_output, em_input) { test_env$population$nfleets <- sum(om_input$fleet_num, om_input$survey_num) test_env$population$nseasons <- 1 test_env$population$nyears <- om_input$nyr + # following line is related to issue #521, may be modified in the future + # https://github.com/NOAA-FIMS/FIMS/issues/521 test_env$population$prop_female <- om_input$proportion.female[1] test_env$population$SetMaturity(test_env$maturity$get_id()) test_env$population$SetGrowth(test_env$ewaa_growth$get_id()) test_env$population$SetRecruitment(test_env$recruitment$get_id()) + # end of setup_fims function, returning test_env return(test_env) } test_that("deterministic test of fims", { + # run function defined above to set up the test environment deterministic_env <- setup_fims( om_input = om_input, om_output = om_output, em_input = em_input ) - # Set-up TMB + # Set-up deterministic_env$fims$CreateTMBModel() - # Create parameter list from Rcpp modules + # CreateTMBModel calls a function in information that loops + # over all the populations and fleets and sets all the pointers parameters <- list(p = deterministic_env$fims$get_fixed()) - par_list <- 1:length(parameters[[1]]) - par_list[2:length(par_list)] <- NA - map <- list(p = factor(par_list)) - - obj <- MakeADFun(data = list(), parameters, DLL = "FIMS", map = map) + # get_fixed function is an Rcpp function that loops over all Rcpp + # modules and returned a vector of parameters being estimated + + # Set up TMB's computational graph + obj <- MakeADFun(data = list(), parameters, DLL = "FIMS") + # Calculate standard errors sdr <- TMB::sdreport(obj) sdr_fixed <- summary(sdr, "fixed") @@ -203,41 +205,46 @@ test_that("deterministic test of fims", { # obj$report() requires parameter list to avoid errors report <- obj$report(obj$par) - # log(R0) + # Compare log(R0) to true value fims_logR0 <- sdr_fixed[1, "Estimate"] expect_gt(fims_logR0, 0.0) expect_equal(fims_logR0, log(om_input$R0)) - # Numbers at age + # Compare numbers at age to true value for (i in 1:length(c(t(om_output$N.age)))) { expect_equal(report$naa[[1]][i], c(t(om_output$N.age))[i]) } - # Biomass + # Compare biomass to true value for (i in 1:length(om_output$biomass.mt)) { expect_equal(report$biomass[[1]][i], om_output$biomass.mt[i]) } - # Spawning biomass + # Compare spawning biomass to true value for (i in 1:length(om_output$SSB)) { expect_equal(report$ssb[[1]][i], om_output$SSB[i]) } - # Recruitment + # Compare recruitment to true value fims_naa <- matrix(report$naa[[1]][1:(om_input$nyr * om_input$nages)], nrow = om_input$nyr, byrow = TRUE ) - for (i in 1:length(om_output$N.age[, 1])) { + # loop over years to compare recruitment by year + for (i in 1:om_input$nyr) { expect_equal(fims_naa[i, 1], om_output$N.age[i, 1]) } + # confirm that recruitment matches the numbers in the first age + # by comparing to fims_naa (what's reported from FIMS) expect_equal( fims_naa[1:om_input$nyr, 1], report$recruitment[[1]][1:om_input$nyr] ) - for (i in 1:length(om_output$N.age[, 1])) { + # confirm that recruitment matches the numbers in the first age + # by comparing to the true values from the OM + for (i in 1:om_input$nyr) { expect_equal(report$recruitment[[1]][i], om_output$N.age[i, 1]) } @@ -250,7 +257,6 @@ test_that("deterministic test of fims", { # Expected catch fims_index <- report$exp_index - # Expect small relative error for deterministic test for (i in 1:length(om_output$L.mt$fleet1)) { expect_equal(fims_index[[1]][i], om_output$L.mt$fleet1[i]) } @@ -264,12 +270,13 @@ test_that("deterministic test of fims", { # Expect 95% of relative error to be within 2*cv expect_lte(sum(fims_object_are > om_input$cv.L$fleet1 * 2.0), length(em_input$L.obs$fleet1) * 0.05) - # Expected catch number at age + # Compare expected catch number at age to true values for (i in 1:length(c(t(om_output$L.age$fleet1)))) { expect_equal(report$cnaa[[1]][i], c(t(om_output$L.age$fleet1))[i]) } # Expected catch number at age in proportion + # QUESTION: Isn't this redundant with the non-proportion test above? fims_cnaa <- matrix(report$cnaa[[1]][1:(om_input$nyr * om_input$nages)], nrow = om_input$nyr, byrow = TRUE ) @@ -280,8 +287,9 @@ test_that("deterministic test of fims", { expect_equal(c(t(fims_cnaa_proportion))[i], c(t(om_cnaa_proportion))[i]) } - # Expected survey index - cwaa <- matrix(report$cwaa[[2]][1:(om_input$nyr * om_input$nages)], nrow = om_input$nyr, byrow = T) + # Expected survey index. + # Using [[2]] because the survey is the 2nd fleet. + cwaa <- matrix(report$cwaa[[2]][1:(om_input$nyr * om_input$nages)], nrow = om_input$nyr, byrow = TRUE) expect_equal(fims_index[[2]], apply(cwaa, 1, sum) * om_output$survey_q$survey1) for (i in 1:length(om_output$survey_index_biomass$survey1)) { @@ -310,7 +318,7 @@ test_that("deterministic test of fims", { for (i in 1:length(c(t(om_cnaa_proportion)))) { expect_equal(c(t(fims_cnaa_proportion))[i], c(t(om_cnaa_proportion))[i]) } - + # clear memory deterministic_env$fims$clear() }) @@ -458,8 +466,8 @@ test_that("estimation test of fims", { ) sdr_naa1_vec <- sdr_report[which(rownames(sdr_report) == "NAA"), 2] sdr_naa1 <- sdr_naa1_vec[seq(1, om_input$nyr * om_input$nages, by = om_input$nages)] - fims_naa1_are <- rep(0, length(om_output$N.age[, 1])) - for (i in 1:length(om_output$N.age[, 1])) { + fims_naa1_are <- rep(0, om_input$nyr) + for (i in 1:om_input$nyr) { fims_naa1_are[i] <- abs(fims_naa[i, 1] - om_output$N.age[i, 1]) # / # om_output$N.age[i, 1] # expect_lte(fims_naa1_are[i], 0.25)