Skip to content

Commit

Permalink
update integration tests with latest code from the main branch
Browse files Browse the repository at this point in the history
  • Loading branch information
@Bai-Li-NOAA @Andrea-Havron-NOAA
Bai-Li-NOAA authored and Andrea-Havron-NOAA committed Aug 27, 2024
1 parent dde3e0e commit 54d97ea
Show file tree
Hide file tree
Showing 3 changed files with 38 additions and 51 deletions.
2 changes: 1 addition & 1 deletion tests/testthat/fixtures/simulate-integration-test-data.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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")
)
)
37 changes: 22 additions & 15 deletions tests/testthat/helper-integration-tests-setup.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,9 @@ setup_and_run_FIMS <- function(iter_id,
em_input_list,
estimation_mode = TRUE,
map = list()) {
# Load operating model data for the current iteration
library(FIMS)

# 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]]
Expand Down Expand Up @@ -79,33 +81,39 @@ 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
# logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is
# 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
Expand Down Expand Up @@ -146,9 +154,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)
Expand Down Expand Up @@ -176,19 +184,17 @@ 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
survey_fleet$SetSelectivity(survey_fleet_selectivity$get_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)
Expand Down Expand Up @@ -251,7 +257,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
))
}

Expand Down
50 changes: 15 additions & 35 deletions tests/testthat/test-integration-fims-estimation.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"]
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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(
Expand All @@ -152,40 +152,19 @@ 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)
fims_logR0 <- sdr_fixed[1, "Estimate"]
# 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(
Expand All @@ -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
))
Expand All @@ -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)
Expand All @@ -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]],
Expand Down

0 comments on commit 54d97ea

Please sign in to comment.