diff --git a/tests/testthat/fixtures/simulate-integration-test-data.R b/tests/testthat/fixtures/simulate-integration-test-data.R index 719822f0..7c0633e4 100644 --- a/tests/testthat/fixtures/simulate-integration-test-data.R +++ b/tests/testthat/fixtures/simulate-integration-test-data.R @@ -45,4 +45,4 @@ for (i in 1:sim_num) { save(om_input_list, om_output_list, em_input_list, file = test_path("fixtures", "integration_test_data.RData") -) +) \ No newline at end of file diff --git a/tests/testthat/helper-integration-tests-setup.R b/tests/testthat/helper-integration-tests-setup.R index c841317e..10176b52 100644 --- a/tests/testthat/helper-integration-tests-setup.R +++ b/tests/testthat/helper-integration-tests-setup.R @@ -48,7 +48,8 @@ setup_and_run_FIMS <- function(iter_id, em_input_list, estimation_mode = TRUE, map = list()) { - # Load operating model data for the current iteration + + # Load operating model data om_input <- om_input_list[[iter_id]] om_output <- om_output_list[[iter_id]] em_input <- em_input_list[[iter_id]] @@ -79,7 +80,7 @@ setup_and_run_FIMS <- function(iter_id, # recruit deviations should enter the model in normal space. # The log is taken in the likelihood calculations # alternative setting: recruitment$log_devs <- rep(0, length(om_input$logR.resid)) - recruitment$log_devs <- methods::new(ParameterVector, om_input$logR.resid[-1], length(om_input$logR.resid[-1]) ) + recruitment$log_devs <- methods::new(ParameterVector, om_input$logR.resid[-1], om_input$nyr-1) recruitment_distribution <- new(TMBDnormDistribution) # set up logR_sd using the normal log_sd parameter @@ -87,25 +88,31 @@ setup_and_run_FIMS <- function(iter_id, # taken before the likelihood calculation recruitment_distribution$log_sd <- new(ParameterVector, 1) recruitment_distribution$log_sd[1]$value <- log(om_input$logR_sd) - recruitment_distribution$log_sd[1]$estimated = FALSE + recruitment_distribution$log_sd[1]$estimated <- FALSE + recruitment_distribution$x <- new(ParameterVector, om_input$nyr) + recruitment_distribution$expected_values <- new(ParameterVector, om_input$nyr) + for (i in 1:om_input$nyr) { + recruitment_distribution$x[i]$value <- 0 + recruitment_distribution$expected_values[i]$value <- 0 + } recruitment_distribution$set_distribution_links("random_effects", recruitment$log_devs$get_id()) - recruitment$estimate_log_devs = TRUE + recruitment$estimate_log_devs <- TRUE # Data 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 - fishing_fleet_index <- new(Index, length(catch)) + fishing_fleet_index <- new(Index, om_input$nyr) fishing_fleet_index$index_data <- catch # set fishing fleet age comp data, need to set dimensions of age comps - fishing_fleet_age_comp <- new(AgeComp, length(catch), om_input$nages) + fishing_fleet_age_comp <- new(AgeComp, om_input$nyr, om_input$nages) fishing_fleet_age_comp$age_comp_data <- c(t(em_input$L.age.obs$fleet1)) * em_input$n.L$fleet1 # repeat for surveys survey_index <- em_input$surveyB.obs$survey1 - survey_fleet_index <- new(Index, length(survey_index)) + survey_fleet_index <- new(Index, om_input$nyr) survey_fleet_index$index_data <- survey_index - survey_fleet_age_comp <- new(AgeComp, length(survey_index), om_input$nages) + survey_fleet_age_comp <- new(AgeComp, om_input$nyr, om_input$nages) survey_fleet_age_comp$age_comp_data <- c(t(em_input$survey.age.obs$survey1)) * em_input$n.survey$survey1 # Growth @@ -146,9 +153,9 @@ setup_and_run_FIMS <- function(iter_id, # Set up fishery index data using the lognormal fishing_fleet_index_distribution <- methods::new(TMBDlnormDistribution) - #lognormal observation error transformed on the log scale + # lognormal observation error transformed on the log scale fishing_fleet_index_distribution$log_logsd <- new(ParameterVector, om_input$nyr) - for(y in 1:om_input$nyr){ + for (y in 1:om_input$nyr) { fishing_fleet_index_distribution$log_logsd[y]$value <- log(sqrt(log(em_input$cv.L$fleet1^2 + 1))) } fishing_fleet_index_distribution$log_logsd$set_all_estimable(FALSE) @@ -176,8 +183,6 @@ setup_and_run_FIMS <- function(iter_id, survey_fleet$is_survey <- TRUE survey_fleet$nages <- om_input$nages survey_fleet$nyears <- om_input$nyr - #survey_fleet$estimate_F <- FALSE - #survey_fleet$random_F <- FALSE survey_fleet$log_q <- log(om_output$survey_q$survey1) survey_fleet$estimate_q <- TRUE survey_fleet$random_q <- FALSE @@ -185,10 +190,10 @@ setup_and_run_FIMS <- function(iter_id, # Set up survey index data using the lognormal survey_fleet_index_distribution <- methods::new(TMBDlnormDistribution) - #lognormal observation error transformed on the log scale + # lognormal observation error transformed on the log scale # sd = sqrt(log(cv^2 + 1)), sd is log transformed survey_fleet_index_distribution$log_logsd <- new(ParameterVector, om_input$nyr) - for(y in 1:om_input$nyr){ + for (y in 1:om_input$nyr) { survey_fleet_index_distribution$log_logsd[y]$value <- log(sqrt(log(em_input$cv.survey$survey1^2 + 1))) } survey_fleet_index_distribution$log_logsd$set_all_estimable(FALSE) @@ -251,7 +256,8 @@ setup_and_run_FIMS <- function(iter_id, opt = opt, report = report, sdr_report = sdr_report, - sdr_fixed = sdr_fixed + sdr_fixed = sdr_fixed, + sdr = sdr )) } diff --git a/tests/testthat/test-integration-fims-estimation.R b/tests/testthat/test-integration-fims-estimation.R index de84a3de..f1e81522 100644 --- a/tests/testthat/test-integration-fims-estimation.R +++ b/tests/testthat/test-integration-fims-estimation.R @@ -14,12 +14,12 @@ test_that("deterministic test of fims", { obj <- result$obj # Calculate standard errors - sdr <- TMB::sdreport(obj) + sdr <- result$sdr sdr_fixed <- result$sdr_fixed -# # Call report using deterministic parameter values -# # obj$report() requires parameter list to avoid errors -# report <- obj$report(obj$par) + # Call report using deterministic parameter values + # obj$report() requires parameter list to avoid errors + report <- result$report # Compare log(R0) to true value fims_logR0 <- sdr_fixed[1, "Estimate"] @@ -99,9 +99,9 @@ test_that("deterministic test of fims", { fims_cnaa_proportion <- fims_cnaa / rowSums(fims_cnaa) om_cnaa_proportion <- om_output_list[[iter_id]]$L.age$fleet1 / rowSums(om_output_list[[iter_id]]$L.age$fleet1) -# for (i in 1:length(c(t(om_cnaa_proportion)))) { -# expect_equal(c(t(fims_cnaa_proportion))[i], c(t(om_cnaa_proportion))[i]) -# } + for (i in 1:length(c(t(om_cnaa_proportion)))) { + expect_equal(c(t(fims_cnaa_proportion))[i], c(t(om_cnaa_proportion))[i]) + } # Expected survey index. # Using [[2]] because the survey is the 2nd fleet. @@ -141,7 +141,7 @@ test_that("deterministic test of fims", { } }) -test_that("nll test of fims", { +test_that("nll test of fims", { iter_id <- 1 result <- setup_and_run_FIMS( @@ -152,26 +152,12 @@ test_that("nll test of fims", { estimation_mode = FALSE ) - parameters <- result$parameters - par_list <- 1:length(parameters[[1]]) - par_list[2:length(par_list)] <- NA - map <- list(p = factor(par_list)) - - result <- setup_and_run_FIMS( - iter_id = iter_id, - om_input_list = om_input_list, - om_output_list = om_output_list, - em_input_list = em_input_list, - estimation_mode = FALSE, - map = map - ) - # Set up TMB's computational graph obj <- result$obj report <- result$report # Calculate standard errors - sdr <- TMB::sdreport(obj) + sdr <- result$sdr sdr_fixed <- result$sdr_fixed # log(R0) @@ -179,13 +165,6 @@ test_that("nll test of fims", { # expect_lte(abs(fims_logR0 - log(om_input$R0)) / log(om_input$R0), 0.0001) expect_equal(fims_logR0, log(om_input_list[[iter_id]]$R0)) - # Call report using deterministic parameter values - # obj$report() requires parameter list to avoid errors - report <- obj$report(obj$par) - obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS", map = map) - jnll <- obj$fn() - - # recruitment likelihood # log_devs is of length nyr-1 rec_nll <- -sum(dnorm( @@ -194,13 +173,13 @@ test_that("nll test of fims", { )) # catch and survey index expected likelihoods - index_nll_fleet <- -sum(dnorm( - log(em_input_list[[iter_id]]$L.obs$fleet1), + index_nll_fleet <- -sum(dlnorm( + em_input_list[[iter_id]]$L.obs$fleet1, log(om_output_list[[iter_id]]$L.mt$fleet1), sqrt(log(em_input_list[[iter_id]]$cv.L$fleet1^2 + 1)), TRUE )) - index_nll_survey <- -sum(dnorm( - log(em_input_list[[iter_id]]$surveyB.obs$survey1), + index_nll_survey <- -sum(dlnorm( + em_input_list[[iter_id]]$surveyB.obs$survey1, log(om_output_list[[iter_id]]$survey_index_biomass$survey1), sqrt(log(em_input_list[[iter_id]]$cv.survey$survey1^2 + 1)), TRUE )) @@ -226,6 +205,7 @@ test_that("nll test of fims", { } age_comp_nll <- age_comp_nll_fleet + age_comp_nll_survey expected_jnll <- rec_nll + index_nll + age_comp_nll + jnll <- report$jnll expect_equal(report$nll_components[1], rec_nll) expect_equal(report$nll_components[2], index_nll_fleet) @@ -249,7 +229,7 @@ test_that("estimation test of fims", { # Compare FIMS results with model comparison project OM values validate_fims( report = result$report, - sdr = TMB::sdreport(result$obj), + sdr = result$sdr, sdr_report = result$sdr_report, om_input = om_input_list[[iter_id]], om_output = om_output_list[[iter_id]],