diff --git a/R/fimsfit.R b/R/fimsfit.R index 8f894031..4e709029 100644 --- a/R/fimsfit.R +++ b/R/fimsfit.R @@ -108,9 +108,11 @@ NULL #' is the returned object from [create_default_parameters()]. #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims methods::setGeneric("get_input", function(x) standardGeneric("get_input")) #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims methods::setMethod("get_input", "FIMSFit", function(x) x@input) #' @return @@ -118,18 +120,22 @@ methods::setMethod("get_input", "FIMSFit", function(x) x@input) #' reportable in the C++ code is returned. #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setGeneric("get_report", function(x) standardGeneric("get_report")) #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setMethod("get_report", "FIMSFit", function(x) x@report) #' @return #' [get_obj()] returns the output from [TMB::MakeADFun()]. #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setGeneric("get_obj", function(x) standardGeneric("get_obj")) #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setMethod("get_obj", "FIMSFit", function(x) x@obj) #' @return @@ -137,9 +143,11 @@ setMethod("get_obj", "FIMSFit", function(x) x@obj) #' in [fit_fims()]. #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setGeneric("get_opt", function(x) standardGeneric("get_opt")) #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setMethod("get_opt", "FIMSFit", function(x) x@opt) #' @return @@ -147,9 +155,11 @@ setMethod("get_opt", "FIMSFit", function(x) x@opt) #' model. #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setGeneric("get_max_gradient", function(x) standardGeneric("get_max_gradient")) #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setMethod("get_max_gradient", "FIMSFit", function(x) x@max_gradient) @@ -157,9 +167,11 @@ setMethod("get_max_gradient", "FIMSFit", function(x) x@max_gradient) #' [get_sdreport()] returns the list from [TMB::sdreport()]. #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setGeneric("get_sdreport", function(x) standardGeneric("get_sdreport")) #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setMethod("get_sdreport", "FIMSFit", function(x) x@sdreport) #' @return @@ -167,9 +179,11 @@ setMethod("get_sdreport", "FIMSFit", function(x) x@sdreport) #' uncertainties from a fitted model. #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setGeneric("get_estimates", function(x) standardGeneric("get_estimates")) #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setMethod("get_estimates", "FIMSFit", function(x) x@estimates) #' @return @@ -178,12 +192,14 @@ setMethod("get_estimates", "FIMSFit", function(x) x@estimates) #' in the model. #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setGeneric( "get_number_of_parameters", function(x) standardGeneric("get_number_of_parameters") ) #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setMethod( "get_number_of_parameters", "FIMSFit", @@ -195,9 +211,11 @@ setMethod( #' seconds as a `difftime` object. #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setGeneric("get_timing", function(x) standardGeneric("get_timing")) #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setMethod("get_timing", "FIMSFit", function(x) x@timing) #' @return @@ -205,9 +223,11 @@ setMethod("get_timing", "FIMSFit", function(x) x@timing) #' the model. #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setGeneric("get_version", function(x) standardGeneric("get_version")) #' @export #' @rdname get_FIMSFit +#' @keywords fit_fims setMethod("get_version", "FIMSFit", function(x) x@version) # methods::setValidity ---- diff --git a/R/fimsframe.R b/R/fimsframe.R index 0b4502ad..6c1666a6 100644 --- a/R/fimsframe.R +++ b/R/fimsframe.R @@ -627,12 +627,13 @@ FIMSFrame <- function(data) { # Get the range of ages displayed in the data to use to specify population # simulation range if ("age" %in% colnames(data)) { - ages <- min(data[["age"]], na.rm = TRUE):max(data[["age"]], na.rm = TRUE) + ages <- sort(unique(data[["age"]])) + ages <- ages[!is.na(ages)] } else { ages <- numeric() } n_ages <- length(ages) - + if ("length" %in% colnames(data)) { lengths <- sort(unique(data[["length"]])) lengths <- lengths[!is.na(lengths)] diff --git a/R/initialize_modules.R b/R/initialize_modules.R index 3934bb5e..323fec43 100644 --- a/R/initialize_modules.R +++ b/R/initialize_modules.R @@ -104,12 +104,17 @@ initialize_module <- function(parameters, data, module_name) { # Set the estimation information for the entire parameter vector module[["age_length_conversion_matrix"]]$set_all_estimable(FALSE) + + module[["age_length_conversion_matrix"]]$set_all_random(FALSE) + } else { + module_fields <- setdiff(module_fields, c( + # Right now we can also remove nlengths because the default is 0 + "nlengths" + )) } module_fields <- setdiff(module_fields, c( "age_length_conversion_matrix", - # Right now we can also remove nlengths because the default is 0 - "nlengths", "proportion_catch_numbers_at_length" )) } @@ -403,12 +408,6 @@ initialize_selectivity <- function(parameters, data, fleet_name) { #' The initialized fleet module as an object. #' @noRd initialize_fleet <- function(parameters, data, fleet_name, linked_ids) { - if (any(is.na(linked_ids[c("selectivity", "index", "age_comp")]))) { - cli::cli_abort(c( - "{.var linked_ids} for {fleet_name} must include 'selectivity', 'index', - and 'age_comp' IDs." - )) - } module <- initialize_module( parameters = parameters, @@ -418,12 +417,23 @@ initialize_fleet <- function(parameters, data, fleet_name, linked_ids) { module$SetSelectivity(linked_ids["selectivity"]) module$SetObservedIndexData(linked_ids["index"]) - module$SetObservedAgeCompData(linked_ids["age_comp"]) fleet_types <- get_data(data) |> dplyr::filter(name == fleet_name) |> dplyr::pull(type) |> unique() + + # Link the observed age composition data to the fleet module using its associated ID + # if the data type includes "age" and if "AgeComp" exists in the data distribution + # specification + if ("age" %in% fleet_types & + "AgeComp" %in% names(parameters[["modules"]][["fleets"]][[fleet_name]][["data_distribution"]])) { + module$SetObservedAgeCompData(linked_ids["age_comp"]) + } + + # Link the observed length composition data to the fleet module using its associated ID + # if the data type includes "length" and if "LengthComp" exists in the data + # distribution specification if ("length" %in% fleet_types & "LengthComp" %in% names(parameters[["modules"]][["fleets"]][[fleet_name]][["data_distribution"]])) { module$SetObservedLengthCompData(linked_ids["length_comp"]) @@ -611,14 +621,8 @@ initialize_fims <- function(parameters, data) { fleet_name = fleet_names[i] ) - fleet_age_comp[[i]] <- initialize_age_comp( - data = data, - fleet_name = fleet_names[i] - ) - fleet_module_ids <- c( index = fleet_index[[i]]$get_id(), - age_comp = fleet_age_comp[[i]]$get_id(), selectivity = fleet_selectivity[[i]]$get_id() ) @@ -627,12 +631,36 @@ initialize_fims <- function(parameters, data) { dplyr::pull(type) |> unique() + # Initialize age composition module if the data type includes "age" and + # if "AgeComp" exists in the data distribution specification + if ("age" %in% fleet_types & + "AgeComp" %in% names(parameters[["modules"]][["fleets"]][[fleet_names[i]]][["data_distribution"]])) { + + # Initialize age composition module for the current fleet + fleet_age_comp[[i]] <- initialize_age_comp( + data = data, + fleet_name = fleet_names[i] + ) + + # Add the module ID for the initialized age composition to the list of fleet module IDs + fleet_module_ids <- c( + fleet_module_ids, + c(age_comp = fleet_age_comp[[i]]$get_id()) + ) + } + + # Initialize length composition module if the data type includes "length" and + # if "LengthComp" exists in the data distribution specification if ("length" %in% fleet_types & "LengthComp" %in% names(parameters[["modules"]][["fleets"]][[fleet_names[i]]][["data_distribution"]])) { + + # Initialize length composition module for the current fleet fleet_length_comp[[i]] <- initialize_length_comp( data = data, fleet_name = fleet_names[i] ) + + # Add the module ID for the initialized length composition to the list of fleet module IDs fleet_module_ids <- c( fleet_module_ids, c(length_comp = fleet_length_comp[[i]]$get_id()) @@ -678,11 +706,14 @@ initialize_fims <- function(parameters, data) { data_type = "index" ) - fleet_agecomp_distribution[[i]] <- initialize_data_distribution( - module = fleet[[i]], - family = multinomial(link = "logit"), - data_type = "agecomp" - ) + if ("age" %in% fleet_types & + "AgeComp" %in% names(parameters[["modules"]][["fleets"]][[fleet_names[i]]][["data_distribution"]])) { + fleet_agecomp_distribution[[i]] <- initialize_data_distribution( + module = fleet[[i]], + family = multinomial(link = "logit"), + data_type = "agecomp" + ) + } if ("length" %in% fleet_types & "LengthComp" %in% names(parameters[["modules"]][["fleets"]][[fleet_names[i]]][["data_distribution"]])) { diff --git a/data-raw/data1.R b/data-raw/data1.R index bbeee22e..9444b606 100644 --- a/data-raw/data1.R +++ b/data-raw/data1.R @@ -62,7 +62,136 @@ ASSAMC::run_om(input_list = sim_input) setwd(working_dir) -load(file.path(main_dir, "sim_data", "output", "OM", paste0("OM", 1, ".RData"))) +# Helper function to calculate length at age using the von Bertalanffy growth model +# a: current age +# Linf: asymptotic average length +# K: Growth coefficient +# a_0: Theoretical age at size zero +AtoL <- function(a, Linf, K, a_0) { + L <- Linf * (1 - exp(-K * (a - a_0))) +} + +# Initialize lists for operating model (OM) and estimation model (EM) inputs and outputs +om_input_list <- om_output_list <- em_input_list <- + vector(mode = "list", length = sim_num) + +# Loop through each simulation to generate length data +for (iter in 1:sim_num) { + # Load the OM data for the current simulation + load(file.path(main_dir, "sim_data", "output", "OM", paste0("OM", iter, ".RData"))) + + # Extract von Bertalanffy growth model parameters from the OM input + Linf <- om_input$Linf + K <- om_input$K + a0 <- om_input$a0 + amax <- max(om_input$ages) + # Define coefficient of variation for length-at-age + cv <- 0.1 + # Extract length-weight coefficient from OM + L2Wa <- om_input$a.lw + # Extract length-weight exponent from OM + L2Wb <- om_input$b.lw + + # Extract age bins from the OM input + ages <- om_input$ages + # Define length bins in intervals of 50 + len_bins <- seq(0, 1100, 50) + + # Create length at age conversion matrix and fill proportions using above + # growth parameters + age_to_length_conversion <- matrix(NA, nrow = length(ages), ncol = length(len_bins)) + for (age in seq_along(ages)) { + # Calculate mean length at age to spread lengths around + mean_length <- AtoL(ages[age], Linf, K, a0) + # mean_length <- AtoLSchnute(ages[age],L1,L2,a1,a2,Ks) + # Calculate the cumulative proportion shorter than each composition length + temp_len_probs <- pnorm(q = len_bins, mean = mean_length, sd = mean_length * cv) + # Reset the first length proportion to zero so the first bin includes all + # density smaller than that bin + temp_len_probs[1] <- 0 + # subtract the offset length probabilities to calculate the proportion in each + # bin. For each length bin the proportion is how many fish are larger than this + # length but shorter than the next bin length. + temp_len_probs <- c(temp_len_probs[-1], 1) - temp_len_probs + age_to_length_conversion[age, ] <- temp_len_probs + } + colnames(age_to_length_conversion) <- len_bins + rownames(age_to_length_conversion) <- ages + + # Loop through each simulation to load the results from the corresponding + # .RData files + # Assign the conversion matrix and other information to the OM input + om_input$lengths <- len_bins + om_input$nlengths <- length(len_bins) + om_input$cv.length_at_age <- cv + om_input$age_to_length_conversion <- age_to_length_conversion + + om_output$L.length <- list() + om_output$survey_length_comp <- list() + om_output$N.length <- matrix(0, nrow = om_input$nyr, ncol = length(len_bins)) + om_output$L.length$fleet1 <- matrix(0, nrow = om_input$nyr, ncol = length(len_bins)) + om_output$survey_length_comp$survey1 <- matrix(0, nrow = om_input$nyr, ncol = length(len_bins)) + + em_input$L.length.obs <- list() + em_input$survey.length.obs <- list() + em_input$L.length.obs$fleet1 <- matrix(0, nrow = om_input$nyr, ncol = length(len_bins)) + em_input$survey.length.obs$survey1 <- matrix(0, nrow = om_input$nyr, ncol = length(len_bins)) + + em_input$lengths <- len_bins + em_input$nlengths <- length(len_bins) + em_input$cv.length_at_age <- cv + em_input$age_to_length_conversion <- age_to_length_conversion + em_input$n.L.lengthcomp$fleet1 <- em_input$n.survey.lengthcomp$survey1 <- 200 + + # Populate length-based outputs for each year, length bin, and age + for (i in seq_along(om_input$year)) { + for (j in seq_along(len_bins)) { + for (k in seq_along(om_input$ages)) { + # Calculate numbers and landings at length for each fleet and survey + om_output$N.length[i, j] <- om_output$N.length[i, j] + + age_to_length_conversion[k, j] * + om_output$N.age[i, k] + + om_output$L.length[[1]][i, j] <- om_output$L.length[[1]][i, j] + + age_to_length_conversion[k, j] * + om_output$L.age[[1]][i, k] + + om_output$survey_length_comp[[1]][i, j] <- om_output$survey_length_comp[[1]][i, j] + + age_to_length_conversion[k, j] * + om_output$survey_age_comp[[1]][i, k] + + em_input$L.length.obs[[1]][i, j] <- em_input$L.length.obs[[1]][i, j] + + age_to_length_conversion[k, j] * + em_input$L.age.obs[[1]][i, k] + + em_input$survey.length.obs[[1]][i, j] <- em_input$survey.length.obs[[1]][i, j] + + age_to_length_conversion[k, j] * + em_input$survey.age.obs[[1]][i, k] + } + } + } + + # Save updated inputs and outputs to file + save( + om_input, om_output, em_input, + file = file.path(main_dir, "sim_data", "output", "OM", paste0("OM", iter, ".RData")) + ) + # Store inputs and outputs in respective lists + om_input_list[[iter]] <- om_input + om_output_list[[iter]] <- om_output + em_input_list[[iter]] <- em_input +} + +# Save all simulations to a single file for {testthat} integration tests +save(om_input_list, om_output_list, em_input_list, + file = testthat::test_path("fixtures", "integration_test_data.RData") +) + +# Load a specific simulation for further processing +sim_id <- 1 +load(file.path(main_dir, "sim_data", "output", "OM", paste0("OM", sim_id, ".RData"))) + +# Return the loaded data returnedom <- list( om_input = om_input, om_output = om_output, @@ -151,7 +280,9 @@ age_data <- rbind( cols = dplyr::starts_with("X"), names_prefix = "X", names_to = "age", - values_to = "value" + values_to = "value", + # Convert the "age" column from strings to integers + names_transform = list(age = as.integer) ) ############################################################################### @@ -180,149 +311,56 @@ weightatage_data <- merge(timingfishery, weightsfishery) ############################################################################### # {FIMS} data ############################################################################### -data1 <- type.convert( - rbind(landings_data, index_data, age_data, weightatage_data), - as.is = TRUE -) +data1 <-rbind(landings_data, index_data, age_data, weightatage_data) # Add new column for length values and set to NA for all milestone 1 data data1 <- data1[, c(1:3, 3:8)] colnames(data1)[4] <- "length" data1[, 4] <- NA -# Growth function values to create length age conversion matrix from model -# comparison project -Linf <- 800 -K <- 0.18 -a0 <- -1.36 -amax <- 12 -cv <- 0.1 - -L2Wa <- 2.5e-08 -L2Wb <- 3 - -AtoL <- function(a, Linf, K, a_0) { - L <- Linf * (1 - exp(-K * (a - a_0))) -} - -ages <- 1:amax -len_bins <- seq(0, 1100, 50) - -# Create length at age conversion matrix and fill proportions using above -# growth parameters -length_age_conversion <- matrix(NA, nrow = length(ages), ncol = length(len_bins)) -for (i in seq_along(ages)) { - # Calculate mean length at age to spread lengths around - mean_length <- AtoL(ages[i], Linf, K, a0) - # mean_length <- AtoLSchnute(ages[i],L1,L2,a1,a2,Ks) - # Calculate the cumulative proportion shorter than each composition length - temp_len_probs <- pnorm(q = len_bins, mean = mean_length, sd = mean_length * cv) - # Reset the first length proportion to zero so the first bin includes all - # density smaller than that bin - temp_len_probs[1] <- 0 - # subtract the offset length probabilities to calculate the proportion in each - # bin. For each length bin the proportion is how many fish are larger than this - # length but shorter than the next bin length. - temp_len_probs <- c(temp_len_probs[-1], 1) - temp_len_probs - length_age_conversion[i, ] <- temp_len_probs -} -colnames(length_age_conversion) <- len_bins -rownames(length_age_conversion) <- ages - # Extract years and fleets from milestone 1 data -start_date <- unique(data1$datestart[data1$type == "weight-at-age"]) -end_date <- unique(data1$dateend[data1$type == "weight-at-age"]) -observers <- unique(data1$name[data1$type == "age"]) +start_date <- timingfishery[["datestart"]] +end_date <- timingfishery[["dateend"]] +observers <- c("fleet1", "survey1") # Create data frame for new fleet and year specific length at age conversion proportions # These are identical across years and fleets in this default example -# length_age_data <- data1[rep(1, length(ages) * -# length(len_bins) * -# length(start_date) * -# length(observers)), ] length_age_data <- data.frame( - matrix(NA, ncol = ncol(data1), nrow = length(ages) * - length(len_bins) * - length(start_date) * - length(observers)) -) -colnames(length_age_data) <- colnames(data1) -row_index <- 1 -for (k in seq_along(start_date)) { - for (l in seq_along(observers)) { - for (i in seq_along(ages)) { - for (j in seq_along(len_bins)) { - length_age_data[row_index, ] <- c( - "age-to-length-conversion", - observers[l], - ages[i], - len_bins[j], - start_date[k], - end_date[k], - length_age_conversion[i, j], - "proportion", - 200 - ) - row_index <- row_index + 1 - } - } - } -} + type = rep("age-to-length-conversion",length(len_bins)*length(ages)*length(observers)*length(start_date)), + name = rep(sort(rep(observers,length(len_bins)*length(ages))),length(start_date)), + age = rep(sort(rep(ages,length(len_bins))),length(observers)*length(start_date)), + length = rep(len_bins,length(ages)*length(observers)*length(start_date)), + datestart = rep(start_date,each=length(len_bins)*length(ages)*length(observers)), + dateend = rep(end_date,each=length(len_bins)*length(ages)*length(observers)), + value = rep(c(t(returnedom$em_input$age_to_length_conversion)),length(observers)*length(start_date)), + unit = rep("proportion",length(len_bins)*length(ages)*length(observers)*length(start_date)), + uncertainty = rep(c(em_input$n.L.lengthcomp$fleet1, em_input$n.survey.lengthcomp$survey1),length(len_bins)*length(ages)*length(start_date))) -# Create a length compostion data frame that will be filled by transforming +# Create a length composition data frame that will be filled by transforming # the age composition data -# length_comp_data <- data1[rep(1, length(len_bins) * -# length(start_date) * -# length(observers)), ] - length_comp_data <- data.frame( - matrix(NA, ncol = ncol(data1), nrow = length(len_bins) * - length(start_date) * - length(observers)) -) -colnames(length_comp_data) <- colnames(data1) - -age_comp_data <- data1[data1$type == "age", ] - -row_index <- 1 -for (k in seq_along(start_date)) { - year_sub <- age_comp_data[age_comp_data$datestart == start_date[k], ] - for (l in seq_along(observers)) { - obs_sub <- year_sub[year_sub$name == observers[l], ] - for (j in seq_along(len_bins)) { - temp_len_props <- length_age_conversion[, j] - temp_len_prob <- sum(temp_len_props * obs_sub$value) - length_comp_data[row_index, ] <- c( - "length", - observers[l], - NA, - len_bins[j], - start_date[k], - end_date[k], - temp_len_prob, - "proportion", - 200 - ) - row_index <- row_index + 1 - } - } -} + type = rep("length",length(len_bins)*length(start_date)*length(observers)), + name = sort(rep(observers,length(len_bins)*length(start_date))), + age = NA, + length = rep(len_bins,length(start_date)*length(observers)), + datestart = rep(start_date,each=length(len_bins)*length(observers)), + dateend = rep(end_date,each=length(len_bins)*length(observers)), + value = c(c(t(returnedom$em_input$L.length.obs$fleet1)), c(t(returnedom$em_input$survey.length.obs$survey1))), + unit = rep("proportion",length(len_bins)*length(start_date)*length(observers)), + uncertainty = rep(c(em_input$n.L.lengthcomp$fleet1, em_input$n.survey.lengthcomp$survey1),length(len_bins)*length(start_date))) # Add the conversion matrix and length composition data to dataframe -# data1 <- rbind(data1,length_comp_data,length_age_data) - -data1 <- type.convert( - rbind(data1, length_comp_data, length_age_data), - as.is = TRUE -) +data1 <- rbind(data1,length_comp_data,length_age_data) write.csv(data1, file.path("FIMS_input_data.csv"), row.names = FALSE ) - # check csv can be read into R well test_read <- read.csv(file.path("FIMS_input_data.csv")) +# TODO: check if we need to do following conversion before running expect_equal() +test_read$datestart <- as.Date(test_read$datestart) +test_read$dateend <- as.Date(test_read$dateend) testthat::expect_equal(test_read, data1) unlink("FIMS_input_data.csv") diff --git a/data/data1.rda b/data/data1.rda index d2591b8d..1ff8968f 100644 Binary files a/data/data1.rda and b/data/data1.rda differ diff --git a/inst/WORDLIST b/inst/WORDLIST index edae1a1b..a1ddb907 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -1,7 +1,11 @@ Arith Beverton +BevertonHoltRecruitment CMD Codecov +devs +DnormDistribution +EWAAgrowth FIMSFit FIMSFrame Github diff --git a/man/get_FIMSFit.Rd b/man/get_FIMSFit.Rd index 4b0e5364..5404d8dd 100644 --- a/man/get_FIMSFit.Rd +++ b/man/get_FIMSFit.Rd @@ -109,3 +109,4 @@ to access objects stored in the available slots. \item \code{\link[=create_default_parameters]{create_default_parameters()}} } } +\keyword{fit_fims} diff --git a/man/m_.Rd b/man/m_.Rd index 85a4c3da..1c8c950a 100644 --- a/man/m_.Rd +++ b/man/m_.Rd @@ -22,11 +22,11 @@ \alias{m_age_to_length_conversion,data.frame-method} \title{Get a vector of data to be passed to a FIMS module from a FIMSFrame object} \usage{ -m_landings(x) +m_landings(x, fleet_name) -\S4method{m_landings}{FIMSFrame}(x) +\S4method{m_landings}{FIMSFrame}(x, fleet_name) -\S4method{m_landings}{data.frame}(x) +\S4method{m_landings}{data.frame}(x, fleet_name) m_index(x, fleet_name) diff --git a/pkgdown/_pkgdown.yml b/pkgdown/_pkgdown.yml index 7f4ce89d..6d1b6d98 100644 --- a/pkgdown/_pkgdown.yml +++ b/pkgdown/_pkgdown.yml @@ -39,6 +39,7 @@ reference: desc: Primary functions used when setting up or running a FIMS model. contents: - starts_with("FIMSFrame") + - starts_with("get_FIMSFrame") - starts_with("m_") - starts_with("create_default_") - update_parameters diff --git a/tests/testthat/fixtures/integration_test_data.RData b/tests/testthat/fixtures/integration_test_data.RData index 8361c1d1..60f924e6 100644 Binary files a/tests/testthat/fixtures/integration_test_data.RData and b/tests/testthat/fixtures/integration_test_data.RData differ diff --git a/tests/testthat/fixtures/simulate-integration-test-data.R b/tests/testthat/fixtures/simulate-integration-test-data.R deleted file mode 100644 index 3f65a5e7..00000000 --- a/tests/testthat/fixtures/simulate-integration-test-data.R +++ /dev/null @@ -1,48 +0,0 @@ -# Install the operating model repo from GitHub -remotes::install_github( - repo = "Bai-Li-NOAA/Age_Structured_Stock_Assessment_Model_Comparison" -) - -working_dir <- getwd() - -maindir <- tempdir() - -# Save the initial OM input using ASSAMC package (sigmaR = 0.4) -model_input <- ASSAMC::save_initial_input() - -# Configure the input parameters for the simulation -sim_num <- 100 -sim_input <- ASSAMC::save_initial_input( - base_case = TRUE, - input_list = model_input, - maindir = maindir, - om_sim_num = sim_num, - keep_sim_num = sim_num, - figure_number = 1, - seed_num = 9924, - case_name = "sim_data" -) - -# Run OM and generate om_input, om_output, and em_input -# using function from the model comparison project -ASSAMC::run_om(input_list = sim_input) - -on.exit(unlink(maindir, recursive = TRUE), add = TRUE) - -setwd(working_dir) -on.exit(setwd(working_dir), add = TRUE) - -# Loop through each simulation to load the results from the corresponding -# .RData files and save them into one file -om_input_list <- om_output_list <- em_input_list <- - vector(mode = "list", length = sim_num) -for (i in 1:sim_num) { - load(file.path(maindir, "sim_data", "output", "OM", paste0("OM", i, ".RData"))) - om_input_list[[i]] <- om_input - om_output_list[[i]] <- om_output - em_input_list[[i]] <- em_input -} - -save(om_input_list, om_output_list, em_input_list, - file = test_path("fixtures", "integration_test_data.RData") -) diff --git a/tests/testthat/helper-integration-tests-setup.R b/tests/testthat/helper-integration-tests-setup.R index cc653875..1d695c49 100644 --- a/tests/testthat/helper-integration-tests-setup.R +++ b/tests/testthat/helper-integration-tests-setup.R @@ -47,17 +47,16 @@ setup_and_run_FIMS_without_wrappers <- function(iter_id, om_output_list, em_input_list, estimation_mode = TRUE, - map = list(), - data) { + map = list()) { # Load operating model data for the current iteration - om_input <- om_input_list[[iter_id]] - om_output <- om_output_list[[iter_id]] - em_input <- em_input_list[[iter_id]] + om_input <- om_input_list[[iter_id]] # Operating model input for the current iteration + om_output <- om_output_list[[iter_id]] # Operating model output for the current iteration + em_input <- em_input_list[[iter_id]] # Estimation model input for the current iteration # Clear any previous FIMS settings clear() - # Data + # Extract fishing fleet landings data (observed) and initialize index module 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 @@ -66,16 +65,13 @@ setup_and_run_FIMS_without_wrappers <- function(iter_id, # set fishing fleet age comp data, need to set dimensions of age comps # Here the new function initializes the object with length nyr*nages fishing_fleet_age_comp <- new(AgeComp, om_input$nyr, om_input$nages) - # Here we fill in the values for the object with the observed length comps for fleet one - # we multiply these proportions by the sample size for likelihood weighting + # Here we fill in the values for the object with the observed age comps for fleet one + # we multiply these proportions by the sample size for likelihood weighting fishing_fleet_age_comp$age_comp_data <- c(t(em_input$L.age.obs$fleet1)) * em_input$n.L$fleet1 - #TODO # set fishing fleet length comp data, need to set dimensions of length comps - fishery_lengthcomp <- FIMS::m_lengthcomp(data, "fleet1") - fishery_age_length_conversion <- FIMS::m_age_to_length_conversion(data, "fleet1") - fishing_fleet_length_comp <- new(LengthComp, om_input$nyr, data@n_lengths) - fishing_fleet_length_comp$length_comp_data <- fishery_lengthcomp * em_input$n.L$fleet1 + fishing_fleet_length_comp <- new(LengthComp, om_input$nyr, om_input$nlengths) + fishing_fleet_length_comp$length_comp_data <- c(t(em_input$L.length.obs$fleet1)) * em_input$n.L.lengthcomp$fleet1 # Fleet # Create the fishing fleet @@ -89,14 +85,18 @@ setup_and_run_FIMS_without_wrappers <- function(iter_id, fishing_fleet_selectivity$slope[1]$is_random_effect <- FALSE fishing_fleet_selectivity$slope[1]$estimated <- TRUE + # Initialize the fishing fleet module fishing_fleet <- new(Fleet) + # Set number of years fishing_fleet$nyears <- om_input$nyr + # Set number of age classes fishing_fleet$nages <- om_input$nages - #TODO - fishing_fleet$nlengths <- data@n_lengths + # Set number of length bins + fishing_fleet$nlengths <- om_input$nlengths fishing_fleet$log_Fmort$resize(om_input$nyr) for (y in 1:om_input$nyr) { + # Log-transform OM fishing mortality fishing_fleet$log_Fmort[y]$value <- log(om_output$f[y]) } fishing_fleet$log_Fmort$set_all_estimable(TRUE) @@ -106,15 +106,14 @@ setup_and_run_FIMS_without_wrappers <- function(iter_id, fishing_fleet$SetSelectivity(fishing_fleet_selectivity$get_id()) fishing_fleet$SetObservedIndexData(fishing_fleet_index$get_id()) fishing_fleet$SetObservedAgeCompData(fishing_fleet_age_comp$get_id()) - - #TODO fishing_fleet$SetObservedLengthCompData(fishing_fleet_length_comp$get_id()) - + # Set up fishery index data using the lognormal fishing_fleet_index_distribution <- methods::new(DlnormDistribution) # lognormal observation error transformed on the log scale fishing_fleet_index_distribution$log_sd$resize(om_input$nyr) for (y in 1:om_input$nyr) { + # Compute lognormal SD from OM coefficient of variation (CV) fishing_fleet_index_distribution$log_sd[y]$value <- log(sqrt(log(em_input$cv.L$fleet1^2 + 1))) } fishing_fleet_index_distribution$log_sd$set_all_estimable(FALSE) @@ -127,32 +126,34 @@ setup_and_run_FIMS_without_wrappers <- function(iter_id, fishing_fleet_agecomp_distribution$set_observed_data(fishing_fleet$GetObservedAgeCompDataID()) fishing_fleet_agecomp_distribution$set_distribution_links("data", fishing_fleet$proportion_catch_numbers_at_age$get_id()) - #TODO change age to length # Set up fishery length composition data using the multinomial fishing_fleet_lengthcomp_distribution <- methods::new(DmultinomDistribution) fishing_fleet_lengthcomp_distribution$set_observed_data(fishing_fleet$GetObservedLengthCompDataID()) fishing_fleet_lengthcomp_distribution$set_distribution_links("data", fishing_fleet$proportion_catch_numbers_at_length$get_id()) - #TODO age to length conversion - #Set age-to-length conversion matrix - fishing_fleet$age_length_conversion_matrix <- new(ParameterVector,fishery_age_length_conversion,om_input$nages*data@n_lengths) + # Set age-to-length conversion matrix + # TODO: If an age_to_length_conversion matrix is provided, the code below + # still executes. Consider adding a check in the Rcpp interface to ensure + # users provide a vector of inputs. + fishing_fleet$age_length_conversion_matrix <- new( + ParameterVector, + c(t(em_input$age_to_length_conversion)), + om_input$nages * om_input$nlengths + ) # Turn off estimation for length-at-age fishing_fleet$age_length_conversion_matrix$set_all_estimable(FALSE) fishing_fleet$age_length_conversion_matrix$set_all_random(FALSE) - # repeat data for surveys + # Repeat similar setup for the survey fleet (e.g., index, age comp, and length comp) + # This includes initializing logistic selectivity, observed data modules, and distribution links. survey_index <- em_input$surveyB.obs$survey1 survey_fleet_index <- new(Index, om_input$nyr) survey_fleet_index$index_data <- survey_index 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 - - #TODO change age to length - survey_lengthcomp <- m_lengthcomp(data, "survey1") - survey_age_length_conversion <- m_age_to_length_conversion(data, "survey1") - survey_fleet_length_comp <- new(LengthComp, om_input$nyr, data@n_lengths) - survey_fleet_length_comp$length_comp_data <- survey_lengthcomp * em_input$n.survey$survey1 - + survey_lengthcomp <- em_input$survey.length.obs$survey1 + survey_fleet_length_comp <- new(LengthComp, om_input$nyr, om_input$nlengths) + survey_fleet_length_comp$length_comp_data <- c(t(survey_lengthcomp)) * em_input$n.survey.lengthcomp$survey1 # Fleet # Create the survey fleet survey_fleet_selectivity <- new(LogisticSelectivity) @@ -169,18 +170,14 @@ setup_and_run_FIMS_without_wrappers <- function(iter_id, survey_fleet$is_survey <- TRUE survey_fleet$nages <- om_input$nages survey_fleet$nyears <- om_input$nyr - - #TODO - survey_fleet$nlengths <- data@n_lengths - + survey_fleet$nlengths <- om_input$nlengths survey_fleet$log_q[1]$value <- log(om_output$survey_q$survey1) + survey_fleet$log_q[1]$estimated <- TRUE survey_fleet$estimate_q <- TRUE survey_fleet$random_q <- FALSE survey_fleet$SetSelectivity(survey_fleet_selectivity$get_id()) survey_fleet$SetObservedIndexData(survey_fleet_index$get_id()) survey_fleet$SetObservedAgeCompData(survey_fleet_age_comp$get_id()) - - #TODO survey_fleet$SetObservedLengthCompData(survey_fleet_length_comp$get_id()) # Set up survey index data using the lognormal @@ -196,16 +193,19 @@ setup_and_run_FIMS_without_wrappers <- function(iter_id, survey_fleet_index_distribution$set_observed_data(survey_fleet$GetObservedIndexDataID()) survey_fleet_index_distribution$set_distribution_links("data", survey_fleet$log_expected_index$get_id()) - # Age composition data + # Age composition distribution survey_fleet_agecomp_distribution <- methods::new(DmultinomDistribution) survey_fleet_agecomp_distribution$set_observed_data(survey_fleet$GetObservedAgeCompDataID()) survey_fleet_agecomp_distribution$set_distribution_links("data", survey_fleet$proportion_catch_numbers_at_age$get_id()) - #TODO + # Length composition distribution survey_fleet_lengthcomp_distribution <- methods::new(DmultinomDistribution) survey_fleet_lengthcomp_distribution$set_observed_data(survey_fleet$GetObservedLengthCompDataID()) - survey_fleet_lengthcomp_distribution$set_distribution_links("data", survey_fleet$proportion_catch_numbers_at_length$get_id()) #Set age to length conversion matrix - survey_fleet$age_length_conversion_matrix <- new(ParameterVector,survey_age_length_conversion,om_input$nages*data@n_lengths) + survey_fleet_lengthcomp_distribution$set_distribution_links("data", survey_fleet$proportion_catch_numbers_at_length$get_id()) # Set age to length conversion matrix + survey_fleet$age_length_conversion_matrix <- new( + ParameterVector, + c(t(em_input$age_to_length_conversion)), + om_input$nages * om_input$nlengths) # Turn off estimation for length-at-age survey_fleet$age_length_conversion_matrix$set_all_estimable(FALSE) survey_fleet$age_length_conversion_matrix$set_all_random(FALSE) @@ -380,120 +380,16 @@ setup_and_run_FIMS_with_wrappers <- function(iter_id, om_input <- om_input_list[[iter_id]] om_output <- om_output_list[[iter_id]] em_input <- em_input_list[[iter_id]] - returnedom <- list( - om_input = om_input, - om_output = om_output, - em_input = em_input - ) + # returnedom <- list( + # om_input = om_input, + # om_output = om_output, + # em_input = em_input + # ) # Clear any previous FIMS settings clear() - # Set up data - cv_2_sd <- function(x) { - sqrt(log(x^2 + 1)) - } - - landings_data <- data.frame( - type = "landings", - name = names(returnedom[["om_output"]]$L.mt)[1], - age = NA, - datestart = as.Date( - paste(returnedom[["om_input"]]$year, 1, 1, sep = "-"), - format = "%Y-%m-%d" - ), - dateend = as.Date( - paste(returnedom[["om_input"]]$year, 12, 31, sep = "-"), - format = "%Y-%m-%d" - ), - value = returnedom[["em_input"]]$L.obs[[1]], - unit = "mt", # metric tons - uncertainty = cv_2_sd(returnedom[["em_input"]]$cv.L[[1]]) - ) - - index_data <- data.frame( - type = "index", - name = names(returnedom[["om_output"]]$survey_index)[1], - age = NA, # Not by age in this case, but there is a by age option. - datestart = as.Date( - paste(returnedom[["om_input"]]$year, 1, 1, sep = "-"), - format = "%Y-%m-%d" - ), - dateend = as.Date( - paste(returnedom[["om_input"]]$year, 1, 1, sep = "-"), - format = "%Y-%m-%d" - ), - value = returnedom[["em_input"]]$surveyB.obs[[1]], - unit = "mt", - uncertainty = cv_2_sd(returnedom[["em_input"]]$cv.survey[[1]]) - ) - - age_data <- rbind( - data.frame( - name = names(returnedom[["em_input"]]$n.L), - returnedom[["em_input"]]$L.age.obs$fleet1, - unit = "proportion", - uncertainty = returnedom[["em_input"]]$n.L$fleet1, - datestart = as.Date( - paste(returnedom[["om_input"]][["year"]], 1, 1, sep = "-"), - "%Y-%m-%d" - ), - dateend = as.Date( - paste(returnedom[["om_input"]][["year"]], 12, 31, sep = "-"), - "%Y-%m-%d" - ) - ), - data.frame( - name = names(returnedom[["om_output"]]$survey_age_comp)[1], - returnedom[["em_input"]]$survey.age.obs[[1]], - unit = "number of fish in proportion", - uncertainty = returnedom[["om_input"]][["n.survey"]][["survey1"]], - datestart = as.Date( - paste(returnedom[["om_input"]][["year"]], 1, 1, sep = "-"), - "%Y-%m-%d" - ), - dateend = as.Date( - paste(returnedom[["om_input"]][["year"]], 1, 1, sep = "-"), - "%Y-%m-%d" - ) - ) - ) |> - dplyr::mutate( - type = "age" - ) |> - tidyr::pivot_longer( - cols = dplyr::starts_with("X"), - names_prefix = "X", - names_to = "age", - values_to = "value" - ) - - timingfishery <- data.frame( - datestart = as.Date( - paste(returnedom[["om_input"]][["year"]], 1, 1, sep = "-"), - "%Y-%m-%d" - ), - dateend = as.Date( - paste(returnedom[["om_input"]][["year"]], 12, 31, sep = "-"), - "%Y-%m-%d" - ) - ) - weightsfishery <- data.frame( - type = "weight-at-age", - name = names(returnedom[["em_input"]]$n.L), - age = seq_along(returnedom[["om_input"]][["W.kg"]]), - value = returnedom[["om_input"]][["W.mt"]], - uncertainty = NA, - unit = "mt" - ) - weightatage_data <- merge(timingfishery, weightsfishery) - - data_dataframe <- type.convert( - rbind(landings_data, index_data, age_data, weightatage_data), - as.is = TRUE - ) - - data <- FIMS::FIMSFrame(data_dataframe) + data <- FIMS::FIMSFrame(data1) # Set up default parameters fleets <- list( @@ -501,14 +397,16 @@ setup_and_run_FIMS_with_wrappers <- function(iter_id, selectivity = list(form = "LogisticSelectivity"), data_distribution = c( Index = "DlnormDistribution", - AgeComp = "DmultinomDistribution" + AgeComp = "DmultinomDistribution", + LengthComp = "DmultinomDistribution" ) ), survey1 = list( selectivity = list(form = "LogisticSelectivity"), data_distribution = c( Index = "DlnormDistribution", - AgeComp = "DmultinomDistribution" + AgeComp = "DmultinomDistribution", + LengthComp = "DmultinomDistribution" ) ) ) diff --git a/tests/testthat/test-initialize_modules.R b/tests/testthat/test-initialize_modules.R index 18c704c9..1a9b2b53 100644 --- a/tests/testthat/test-initialize_modules.R +++ b/tests/testthat/test-initialize_modules.R @@ -5,7 +5,8 @@ fleet1 <- survey1 <- list( selectivity = list(form = "LogisticSelectivity"), data_distribution = c( Index = "DlnormDistribution", - AgeComp = "DmultinomDistribution" + AgeComp = "DmultinomDistribution", + LengthComp = "DmultinomDistribution" ) ) diff --git a/tests/testthat/test-integration-fims-estimation-with-wrappers.R b/tests/testthat/test-integration-fims-estimation-with-wrappers.R index 3e784e9a..4cb56210 100644 --- a/tests/testthat/test-integration-fims-estimation-with-wrappers.R +++ b/tests/testthat/test-integration-fims-estimation-with-wrappers.R @@ -1,67 +1,5 @@ load(test_path("fixtures", "integration_test_data.RData")) -fleets <- list( - fleet1 = list( - selectivity = list(form = "LogisticSelectivity"), - data_distribution = c( - Index = "DlnormDistribution", - AgeComp = "DmultinomDistribution" - ) - ), - survey1 = list( - selectivity = list(form = "LogisticSelectivity"), - data_distribution = c( - Index = "DlnormDistribution", - AgeComp = "DmultinomDistribution" - ) - ) -) - -data("data1") -data <- FIMS::FIMSFrame(data1) -default_parameters <- data1 |> - FIMS::FIMSFrame() |> - create_default_parameters( - fleets = fleets, - recruitment = list( - form = "BevertonHoltRecruitment", - process_distribution = c(log_devs = "DnormDistribution") - ), - growth = list(form = "EWAAgrowth"), - maturity = list(form = "LogisticMaturity") - ) - -modified_parameters <- list( - fleet1 = list( - Fleet.log_Fmort.value = log(om_output_list[[1]]$f) - ), - survey1 = list( - LogisticSelectivity.inflection_point.value = 1.5, - LogisticSelectivity.slope.value = 2, - Fleet.log_q.value = log(om_output_list[[1]]$survey_q$survey1) - ), - recruitment = list( - BevertonHoltRecruitment.log_rzero.value = log(om_input_list[[1]]$R0), - BevertonHoltRecruitment.log_devs.value = om_input_list[[1]]$logR.resid[-1], - BevertonHoltRecruitment.log_devs.estimated = FALSE, - DnormDistribution.log_sd.value = om_input_list[[1]]$logR_sd - ), - maturity = list( - LogisticMaturity.inflection_point.value = om_input_list[[1]]$A50.mat, - LogisticMaturity.inflection_point.estimated = FALSE, - LogisticMaturity.slope.value = om_input_list[[1]]$slope.mat, - LogisticMaturity.slope.estimated = FALSE - ), - population = list( - Population.log_init_naa.value = log(om_output_list[[1]]$N.age[1, ]) - ) -) - -parameters <- default_parameters |> - update_parameters( - modified_parameters = modified_parameters - ) - test_that("deterministic test of fims", { iter_id <- 1 @@ -259,15 +197,41 @@ 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 + + # length comp likelihoods + # TODO: the commented-out code below is not working yet + # fishing_lengthcomp_observed <- em_input_list[[iter_id]]$L.length.obs$fleet1 + # fishing_lengthcomp_expected <- om_output_list[[iter_id]]$L.length$fleet1 / rowSums(om_output_list[[iter_id]]$L.length$fleet1) + # survey_lengthcomp_observed <- em_input_list[[iter_id]]$survey.length.obs$survey1 + # survey_lengthcomp_expected <- om_output_list[[iter_id]]$survey_length_comp$survey1 / rowSums(om_output_list[[iter_id]]$survey_length_comp$survey1) + # lengthcomp_nll_fleet <- lengthcomp_nll_survey <- 0 + # for (y in 1:om_input_list[[iter_id]]$nyr) { + # + # lengthcomp_nll_fleet <- lengthcomp_nll_fleet - + # dmultinom( + # fishing_lengthcomp_observed[y, ] * em_input_list[[iter_id]]$n.L.lengthcomp$fleet1, em_input_list[[iter_id]]$n.L.lengthcomp$fleet1, + # fishing_lengthcomp_expected[y, ], TRUE + # ) + # + # lengthcomp_nll_survey <- lengthcomp_nll_survey - + # dmultinom( + # survey_lengthcomp_observed[y, ] * em_input_list[[iter_id]]$n.survey.lengthcomp$survey1, em_input_list[[iter_id]]$n.survey.lengthcomp$survey1, + # survey_lengthcomp_expected[y, ], TRUE + # ) + # } + # lengthcomp_nll <- lengthcomp_nll_fleet + lengthcomp_nll_survey + # + # expected_jnll <- rec_nll + index_nll + age_comp_nll + lengthcomp_nll + # jnll <- report$jnll expect_equal(report$nll_components[1], rec_nll) expect_equal(report$nll_components[2], index_nll_fleet) expect_equal(report$nll_components[3], age_comp_nll_fleet) - expect_equal(report$nll_components[4], index_nll_survey) - expect_equal(report$nll_components[5], age_comp_nll_survey) - expect_equal(jnll, expected_jnll) + # expect_equal(report$nll_components[4], lengthcomp_nll_fleet) + expect_equal(report$nll_components[5], index_nll_survey) + expect_equal(report$nll_components[6], age_comp_nll_survey) + # expect_equal(report$nll_components[7], lengthcomp_nll_survey) + # expect_equal(jnll, expected_jnll) }) test_that("estimation test of fims using wrapper functions", { @@ -293,18 +257,174 @@ test_that("estimation test of fims using wrapper functions", { ) }) -test_that("estimation test of fims using high-level wrappers", { +test_that("estimation test with age and length comp using wrappers",{ + # Load operating model data for the current iteration + iter_id <- 1 + om_input <- om_input_list[[iter_id]] + om_output <- om_output_list[[iter_id]] + em_input <- em_input_list[[iter_id]] + + fims_data <- FIMS::FIMSFrame(data1) + + # Clear any previous FIMS settings + clear() + + fleets <- list( + fleet1 = list( + selectivity = list(form = "LogisticSelectivity"), + data_distribution = c( + Index = "DlnormDistribution", + AgeComp = "DmultinomDistribution", + LengthComp = "DmultinomDistribution" + ) + ), + survey1 = list( + selectivity = list(form = "LogisticSelectivity"), + data_distribution = c( + Index = "DlnormDistribution", + AgeComp = "DmultinomDistribution", + LengthComp = "DmultinomDistribution" + ) + ) + ) + + lengthcomp_parameters <- fims_data |> + create_default_parameters( + fleets = fleets, + recruitment = list( + form = "BevertonHoltRecruitment", + process_distribution = c(log_devs = "DnormDistribution") + ), + growth = list(form = "EWAAgrowth"), + maturity = list(form = "LogisticMaturity") + ) + + modified_parameters <- list( + fleet1 = list( + Fleet.log_Fmort.value = log(om_output_list[[1]]$f) + ), + survey1 = list( + LogisticSelectivity.inflection_point.value = 1.5, + LogisticSelectivity.slope.value = 2, + Fleet.log_q.value = log(om_output_list[[1]]$survey_q$survey1) + ), + recruitment = list( + BevertonHoltRecruitment.log_rzero.value = log(om_input_list[[1]]$R0), + BevertonHoltRecruitment.log_devs.value = om_input_list[[1]]$logR.resid[-1], + BevertonHoltRecruitment.log_devs.estimated = FALSE, + DnormDistribution.log_sd.value = om_input_list[[1]]$logR_sd + ), + maturity = list( + LogisticMaturity.inflection_point.value = om_input_list[[1]]$A50.mat, + LogisticMaturity.inflection_point.estimated = FALSE, + LogisticMaturity.slope.value = om_input_list[[1]]$slope.mat, + LogisticMaturity.slope.estimated = FALSE + ), + population = list( + Population.log_init_naa.value = log(om_output_list[[1]]$N.age[1, ]) + ) + ) + + parameters <- lengthcomp_parameters |> + update_parameters( + modified_parameters = modified_parameters + ) + + parameter_list <- initialize_fims( + parameters = parameters, + data = fims_data + ) + fit <- fit_fims(parameter_list, optimize = TRUE) + + clear() + + validate_fims( + report = fit@report, + sdr = fit@estimates, + sdr_report = fit@estimates, + om_input = om_input_list[[iter_id]], + om_output = om_output_list[[iter_id]], + em_input = em_input_list[[iter_id]], + use_fimsfit = TRUE + ) +}) + +test_that("estimation test with length comp using wrappers",{ # Load operating model data for the current iteration iter_id <- 1 om_input <- om_input_list[[iter_id]] om_output <- om_output_list[[iter_id]] em_input <- em_input_list[[iter_id]] + fims_data <- data1 |> + dplyr::filter(type != "age") |> + FIMS::FIMSFrame() + # Clear any previous FIMS settings clear() + + fleets <- list( + fleet1 = list( + selectivity = list(form = "LogisticSelectivity"), + data_distribution = c( + Index = "DlnormDistribution", + LengthComp = "DmultinomDistribution" + ) + ), + survey1 = list( + selectivity = list(form = "LogisticSelectivity"), + data_distribution = c( + Index = "DlnormDistribution", + LengthComp = "DmultinomDistribution" + ) + ) + ) + + lengthcomp_parameters <- fims_data |> + create_default_parameters( + fleets = fleets, + recruitment = list( + form = "BevertonHoltRecruitment", + process_distribution = c(log_devs = "DnormDistribution") + ), + growth = list(form = "EWAAgrowth"), + maturity = list(form = "LogisticMaturity") + ) + + modified_parameters <- list( + fleet1 = list( + Fleet.log_Fmort.value = log(om_output_list[[1]]$f) + ), + survey1 = list( + LogisticSelectivity.inflection_point.value = 1.5, + LogisticSelectivity.slope.value = 2, + Fleet.log_q.value = log(om_output_list[[1]]$survey_q$survey1) + ), + recruitment = list( + BevertonHoltRecruitment.log_rzero.value = log(om_input_list[[1]]$R0), + BevertonHoltRecruitment.log_devs.value = om_input_list[[1]]$logR.resid[-1], + BevertonHoltRecruitment.log_devs.estimated = FALSE, + DnormDistribution.log_sd.value = om_input_list[[1]]$logR_sd + ), + maturity = list( + LogisticMaturity.inflection_point.value = om_input_list[[1]]$A50.mat, + LogisticMaturity.inflection_point.estimated = FALSE, + LogisticMaturity.slope.value = om_input_list[[1]]$slope.mat, + LogisticMaturity.slope.estimated = FALSE + ), + population = list( + Population.log_init_naa.value = log(om_output_list[[1]]$N.age[1, ]) + ) + ) + + parameters <- lengthcomp_parameters |> + update_parameters( + modified_parameters = modified_parameters + ) + parameter_list <- initialize_fims( parameters = parameters, - data = data + data = fims_data ) fit <- fit_fims(parameter_list, optimize = TRUE) diff --git a/tests/testthat/test-integration-fims-estimation-without-wrappers.R b/tests/testthat/test-integration-fims-estimation-without-wrappers.R index b0935454..1db1cf42 100644 --- a/tests/testthat/test-integration-fims-estimation-without-wrappers.R +++ b/tests/testthat/test-integration-fims-estimation-without-wrappers.R @@ -1,339 +1,334 @@ -load(test_path("fixtures", "integration_test_data.RData")) -data("data1") -data <- FIMS::FIMSFrame(data1) - -test_that("deterministic test of fims", { - iter_id <- 1 - result <- setup_and_run_FIMS_without_wrappers( - 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, - data = data - ) - - # Set up TMB's computational graph - obj <- result$obj - - # Calculate standard errors - sdr <- result$sdr - sdr_fixed <- result$sdr_fixed - - # 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[35, "Estimate"] - expect_gt(fims_logR0, 0.0) - expect_equal(fims_logR0, log(om_input_list[[iter_id]]$R0)) - - # Compare numbers at age to true value - for (i in 1:length(c(t(om_output_list[[iter_id]]$N.age)))) { - expect_equal(report$naa[[1]][i], c(t(om_output_list[[iter_id]]$N.age))[i]) - } - - # Compare biomass to true value - for (i in 1:length(om_output_list[[iter_id]]$biomass.mt)) { - expect_equal(report$biomass[[1]][i], om_output_list[[iter_id]]$biomass.mt[i]) - } - - # Compare spawning biomass to true value - for (i in 1:length(om_output_list[[iter_id]]$SSB)) { - expect_equal(report$ssb[[1]][i], om_output_list[[iter_id]]$SSB[i]) - } - - # Compare recruitment to true value - fims_naa <- matrix(report$naa[[1]][1:(om_input_list[[iter_id]]$nyr * om_input_list[[iter_id]]$nages)], - nrow = om_input_list[[iter_id]]$nyr, byrow = TRUE - ) - - # loop over years to compare recruitment by year - for (i in 1:om_input_list[[iter_id]]$nyr) { - expect_equal(fims_naa[i, 1], om_output_list[[iter_id]]$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_list[[iter_id]]$nyr, 1], - report$recruitment[[1]][1:om_input_list[[iter_id]]$nyr] - ) - - # 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_list[[iter_id]]$nyr) { - expect_equal(report$recruitment[[1]][i], om_output_list[[iter_id]]$N.age[i, 1]) - } - - # recruitment log_devs (fixed at initial "true" values) - # the initial value of om_input$logR.resid is dropped from the model - expect_equal(report$log_recruit_dev[[1]], om_input_list[[iter_id]]$logR.resid[-1]) - - # F (fixed at initial "true" values) - expect_equal(report$F_mort[[1]], om_output_list[[iter_id]]$f) - - # Expected catch - fims_index <- report$exp_index - for (i in 1:length(om_output_list[[iter_id]]$L.mt$fleet1)) { - expect_equal(fims_index[[1]][i], om_output_list[[iter_id]]$L.mt$fleet1[i]) - } - - # Expect small relative error for deterministic test - fims_object_are <- rep(0, length(em_input_list[[iter_id]]$L.obs$fleet1)) - for (i in 1:length(em_input_list[[iter_id]]$L.obs$fleet1)) { - fims_object_are[i] <- abs(fims_index[[1]][i] - em_input_list[[iter_id]]$L.obs$fleet1[i]) / em_input_list[[iter_id]]$L.obs$fleet1[i] - } - - # Expect 95% of relative error to be within 2*cv - expect_lte(sum(fims_object_are > om_input_list[[iter_id]]$cv.L$fleet1 * 2.0), length(em_input_list[[iter_id]]$L.obs$fleet1) * 0.05) - - # Compare expected catch number at age to true values - for (i in 1:length(c(t(om_output_list[[iter_id]]$L.age$fleet1)))) { - expect_equal(report$cnaa[[1]][i], c(t(om_output_list[[iter_id]]$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_list[[iter_id]]$nyr * om_input_list[[iter_id]]$nages)], - nrow = om_input_list[[iter_id]]$nyr, byrow = TRUE - ) - 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]) - } - - # Expected survey index. - # Using [[2]] because the survey is the 2nd fleet. - cwaa <- matrix(report$cwaa[[2]][1:(om_input_list[[iter_id]]$nyr * om_input_list[[iter_id]]$nages)], - nrow = om_input_list[[iter_id]]$nyr, byrow = TRUE - ) - expect_equal(fims_index[[2]], apply(cwaa, 1, sum) * om_output_list[[iter_id]]$survey_q$survey1) - - for (i in 1:length(om_output_list[[iter_id]]$survey_index_biomass$survey1)) { - expect_equal(fims_index[[2]][i], om_output_list[[iter_id]]$survey_index_biomass$survey1[i]) - } - - fims_object_are <- rep(0, length(em_input_list[[iter_id]]$surveyB.obs$survey1)) - for (i in 1:length(em_input_list[[iter_id]]$survey.obs$survey1)) { - fims_object_are[i] <- abs(fims_index[[2]][i] - em_input_list[[iter_id]]$surveyB.obs$survey1[i]) / em_input_list[[iter_id]]$surveyB.obs$survey1[i] - } - # Expect 95% of relative error to be within 2*cv - expect_lte( - sum(fims_object_are > om_input_list[[iter_id]]$cv.survey$survey1 * 2.0), - length(em_input_list[[iter_id]]$surveyB.obs$survey1) * 0.05 - ) - - # Expected catch number at age in proportion - fims_cnaa <- matrix(report$cnaa[[2]][1:(om_input_list[[iter_id]]$nyr * om_input_list[[iter_id]]$nages)], - nrow = om_input_list[[iter_id]]$nyr, byrow = TRUE - ) - - for (i in 1:length(c(t(om_output_list[[iter_id]]$survey_age_comp$survey1)))) { - expect_equal(report$cnaa[[2]][i], c(t(om_output_list[[iter_id]]$survey_age_comp$survey1))[i]) - } - - fims_cnaa_proportion <- fims_cnaa / rowSums(fims_cnaa) - om_cnaa_proportion <- om_output_list[[iter_id]]$survey_age_comp$survey1 / rowSums(om_output_list[[iter_id]]$survey_age_comp$survey1) - - for (i in 1:length(c(t(om_cnaa_proportion)))) { - expect_equal(c(t(fims_cnaa_proportion))[i], c(t(om_cnaa_proportion))[i]) - } -}) - -test_that("nll test of fims", { - iter_id <- 1 - - result <- setup_and_run_FIMS_without_wrappers( - 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, - data = data - ) - - # Set up TMB's computational graph - obj <- result$obj - report <- result$report - - # Calculate standard errors - sdr <- result$sdr - sdr_fixed <- result$sdr_fixed - - # log(R0) - fims_logR0 <- sdr_fixed[35, "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)) - - # recruitment likelihood - # log_devs is of length nyr-1 - rec_nll <- -sum(dnorm( - om_input_list[[iter_id]]$logR.resid[-1], rep(0, om_input_list[[iter_id]]$nyr - 1), - om_input_list[[iter_id]]$logR_sd, TRUE - )) - - # catch and survey index expected likelihoods - 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(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 - )) - index_nll <- index_nll_fleet + index_nll_survey - # age comp likelihoods - fishing_acomp_observed <- em_input_list[[iter_id]]$L.age.obs$fleet1 - fishing_acomp_expected <- om_output_list[[iter_id]]$L.age$fleet1 / rowSums(om_output_list[[iter_id]]$L.age$fleet1) - survey_acomp_observed <- em_input_list[[iter_id]]$survey.age.obs$survey1 - survey_acomp_expected <- om_output_list[[iter_id]]$survey_age_comp$survey1 / rowSums(om_output_list[[iter_id]]$survey_age_comp$survey1) - age_comp_nll_fleet <- age_comp_nll_survey <- 0 - for (y in 1:om_input_list[[iter_id]]$nyr) { - age_comp_nll_fleet <- age_comp_nll_fleet - - dmultinom( - fishing_acomp_observed[y, ] * em_input_list[[iter_id]]$n.L$fleet1, em_input_list[[iter_id]]$n.L$fleet1, - fishing_acomp_expected[y, ], TRUE - ) - - age_comp_nll_survey <- age_comp_nll_survey - - dmultinom( - survey_acomp_observed[y, ] * em_input_list[[iter_id]]$n.survey$survey1, em_input_list[[iter_id]]$n.survey$survey1, - survey_acomp_expected[y, ], TRUE - ) - } - age_comp_nll <- age_comp_nll_fleet + age_comp_nll_survey - - # TODO: length comp likelihoods - # fishing_lcomp_observed <- em_input_list[[iter_id]]$L.length.obs$fleet1 - # fishing_lcomp_expected <- om_output_list[[iter_id]]$L.length$fleet1 / rowSums(om_output_list[[iter_id]]$L.length$fleet1) - # survey_lcomp_observed <- em_input_list[[iter_id]]$survey.length.obs$survey1 - # survey_lcomp_expected <- om_output_list[[iter_id]]$survey_length_comp$survey1 / rowSums(om_output_list[[iter_id]]$survey_length_comp$survey1) - # length_comp_nll_fleet <- length_comp_nll_survey <- 0 - # for (y in 1:om_input_list[[iter_id]]$nyr) { - # age_comp_nll_fleet <- age_comp_nll_fleet - - # dmultinom( - # fishing_lcomp_observed[y, ] * em_input_list[[iter_id]]$n.L$fleet1, em_input_list[[iter_id]]$n.L$fleet1, - # fishing_lcomp_expected[y, ], TRUE - # ) - - # length_comp_nll_survey <- length_comp_nll_survey - - # dmultinom( - # survey_lcomp_observed[y, ] * em_input_list[[iter_id]]$n.survey$survey1, em_input_list[[iter_id]]$n.survey$survey1, - # survey_lcomp_expected[y, ], TRUE - # ) - # } - # 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) - expect_equal(report$nll_components[3], age_comp_nll_fleet) - # expect_equal(report$nll_components[4], length_comp_nll_fleet) - expect_equal(report$nll_components[5], index_nll_survey) - expect_equal(report$nll_components[6], age_comp_nll_survey) - # expect_equal(report$nll_components[7], length_comp_nll_survey) - # expect_equal(jnll, expected_jnll) -}) - -test_that("estimation test of fims", { - # Initialize the iteration identifier and run FIMS with the 1st set of OM values - iter_id <- 1 - result <- setup_and_run_FIMS_without_wrappers( - iter_id = iter_id, - om_input_list = om_input_list, - om_output_list = om_output_list, - em_input_list = em_input_list, - estimation_mode = TRUE, - data = data - ) - - # Compare FIMS results with model comparison project OM values - validate_fims( - report = result$report, - sdr = result$sdr, - sdr_report = result$sdr_report, - om_input = om_input_list[[iter_id]], - om_output = om_output_list[[iter_id]], - em_input = em_input_list[[iter_id]] - ) -}) - -test_that("run FIMS with missing values", { - # Initialize the iteration identifier - iter_id <- 1 - - # Define the NA (missing value) placeholder and the index where it will be inserted - na_value <- -999 - na_index <- 2 - - # Introduce a missing value into the survey observations for the estimation model input - em_input_list[[iter_id]]$surveyB.obs$survey1[na_index] <- na_value - - # Run the FIMS setup and execution function - result <- setup_and_run_FIMS_without_wrappers( - iter_id = iter_id, - om_input_list = om_input_list, - om_output_list = om_output_list, - em_input_list = em_input_list, - estimation_mode = TRUE, - data = data - ) - - # Validate that the result report is not null - expect_false(is.null(result$report)) - - # Obtain the gradient and Hessian matrix - g <- as.numeric(result$obj$gr(result$opt$par)) - h <- optimHess(result$opt$par, fn = result$obj$fn, gr = result$obj$gr) - result$opt$par <- result$opt$par - solve(h, g) - - # Obtain the maximum absolute gradient to check convergence - # Ensure that the maximum gradient is less than or equal to - # the specified tolerance (0.0001) - max_gradient <- max(abs(result$obj$gr(result$opt$par))) - expect_lte(max_gradient, 0.0001) -}) - -test_that("agecomp in proportion works", { - # Initialize the iteration identifier - iter_id <- 1 - - # Store the original values of the number of landings observations and - # survey observations - n.L_original <- om_input_list[[iter_id]]$n.L$fleet1 - n.survey_original <- om_input_list[[iter_id]]$n.survey$survey1 - - # Set the number of landings observations and survey observations to 1 - om_input_list[[iter_id]]$n.L$fleet1 <- 1 - om_input_list[[iter_id]]$n.survey$survey1 <- 1 - on.exit(om_input_list[[iter_id]]$n.L$fleet1 <- n.L_original, add = TRUE) - on.exit(om_input_list[[iter_id]]$n.survey$survey1 <- n.survey_original, add = TRUE) - - # Run the FIMS setup and execution function - result <- setup_and_run_FIMS_without_wrappers( - iter_id = iter_id, - om_input_list = om_input_list, - om_output_list = om_output_list, - em_input_list = em_input_list, - estimation_mode = TRUE, - data = data - ) - - # Compare FIMS results with model comparison project OM values - validate_fims( - report = result$report, - sdr = TMB::sdreport(result$obj), - sdr_report = result$sdr_report, - om_input = om_input_list[[iter_id]], - om_output = om_output_list[[iter_id]], - em_input = em_input_list[[iter_id]] - ) -}) +# Load test data for integration testing +# The test data is stored as an RData file in the tests/testthat/fixtures folder, +# which contains 100 sets of simulated data using {ASSAMC} from +# https://github.com/Bai-Li-NOAA/Age_Structured_Stock_Assessment_Model_Comparison. +load(test_path("fixtures", "integration_test_data.RData")) + +# Initialize the iteration identifier and run FIMS with the 1st set of OM values +iter_id <- 1 + +test_that("deterministic test of fims", { + + result <- setup_and_run_FIMS_without_wrappers( + 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 + ) + + # Set up TMB's computational graph + obj <- result$obj + + # Calculate standard errors + sdr <- result$sdr + sdr_fixed <- result$sdr_fixed + + # 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[36, "Estimate"] + expect_gt(fims_logR0, 0.0) + expect_equal(fims_logR0, log(om_input_list[[iter_id]]$R0)) + + # Compare numbers at age to true value + for (i in 1:length(c(t(om_output_list[[iter_id]]$N.age)))) { + expect_equal(report$naa[[1]][i], c(t(om_output_list[[iter_id]]$N.age))[i]) + } + + # Compare biomass to true value + for (i in 1:length(om_output_list[[iter_id]]$biomass.mt)) { + expect_equal(report$biomass[[1]][i], om_output_list[[iter_id]]$biomass.mt[i]) + } + + # Compare spawning biomass to true value + for (i in 1:length(om_output_list[[iter_id]]$SSB)) { + expect_equal(report$ssb[[1]][i], om_output_list[[iter_id]]$SSB[i]) + } + + # Compare recruitment to true value + fims_naa <- matrix(report$naa[[1]][1:(om_input_list[[iter_id]]$nyr * om_input_list[[iter_id]]$nages)], + nrow = om_input_list[[iter_id]]$nyr, byrow = TRUE + ) + + # loop over years to compare recruitment by year + for (i in 1:om_input_list[[iter_id]]$nyr) { + expect_equal(fims_naa[i, 1], om_output_list[[iter_id]]$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_list[[iter_id]]$nyr, 1], + report$recruitment[[1]][1:om_input_list[[iter_id]]$nyr] + ) + + # 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_list[[iter_id]]$nyr) { + expect_equal(report$recruitment[[1]][i], om_output_list[[iter_id]]$N.age[i, 1]) + } + + # recruitment log_devs (fixed at initial "true" values) + # the initial value of om_input$logR.resid is dropped from the model + expect_equal(report$log_recruit_dev[[1]], om_input_list[[iter_id]]$logR.resid[-1]) + + # F (fixed at initial "true" values) + expect_equal(report$F_mort[[1]], om_output_list[[iter_id]]$f) + + # Expected catch + fims_index <- report$exp_index + for (i in 1:length(om_output_list[[iter_id]]$L.mt$fleet1)) { + expect_equal(fims_index[[1]][i], om_output_list[[iter_id]]$L.mt$fleet1[i]) + } + + # Expect small relative error for deterministic test + fims_object_are <- rep(0, length(em_input_list[[iter_id]]$L.obs$fleet1)) + for (i in 1:length(em_input_list[[iter_id]]$L.obs$fleet1)) { + fims_object_are[i] <- abs(fims_index[[1]][i] - em_input_list[[iter_id]]$L.obs$fleet1[i]) / em_input_list[[iter_id]]$L.obs$fleet1[i] + } + + # Expect 95% of relative error to be within 2*cv + expect_lte(sum(fims_object_are > om_input_list[[iter_id]]$cv.L$fleet1 * 2.0), length(em_input_list[[iter_id]]$L.obs$fleet1) * 0.05) + + # Compare expected catch number at age to true values + for (i in 1:length(c(t(om_output_list[[iter_id]]$L.age$fleet1)))) { + expect_equal(report$cnaa[[1]][i], c(t(om_output_list[[iter_id]]$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_list[[iter_id]]$nyr * om_input_list[[iter_id]]$nages)], + nrow = om_input_list[[iter_id]]$nyr, byrow = TRUE + ) + 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]) + } + + # Expected survey index. + # Using [[2]] because the survey is the 2nd fleet. + cwaa <- matrix(report$cwaa[[2]][1:(om_input_list[[iter_id]]$nyr * om_input_list[[iter_id]]$nages)], + nrow = om_input_list[[iter_id]]$nyr, byrow = TRUE + ) + expect_equal(fims_index[[2]], apply(cwaa, 1, sum) * om_output_list[[iter_id]]$survey_q$survey1) + + for (i in 1:length(om_output_list[[iter_id]]$survey_index_biomass$survey1)) { + expect_equal(fims_index[[2]][i], om_output_list[[iter_id]]$survey_index_biomass$survey1[i]) + } + + fims_object_are <- rep(0, length(em_input_list[[iter_id]]$surveyB.obs$survey1)) + for (i in 1:length(em_input_list[[iter_id]]$survey.obs$survey1)) { + fims_object_are[i] <- abs(fims_index[[2]][i] - em_input_list[[iter_id]]$surveyB.obs$survey1[i]) / em_input_list[[iter_id]]$surveyB.obs$survey1[i] + } + # Expect 95% of relative error to be within 2*cv + expect_lte( + sum(fims_object_are > om_input_list[[iter_id]]$cv.survey$survey1 * 2.0), + length(em_input_list[[iter_id]]$surveyB.obs$survey1) * 0.05 + ) + + # Expected catch number at age in proportion + fims_cnaa <- matrix(report$cnaa[[2]][1:(om_input_list[[iter_id]]$nyr * om_input_list[[iter_id]]$nages)], + nrow = om_input_list[[iter_id]]$nyr, byrow = TRUE + ) + + for (i in 1:length(c(t(om_output_list[[iter_id]]$survey_age_comp$survey1)))) { + expect_equal(report$cnaa[[2]][i], c(t(om_output_list[[iter_id]]$survey_age_comp$survey1))[i]) + } + + fims_cnaa_proportion <- fims_cnaa / rowSums(fims_cnaa) + om_cnaa_proportion <- om_output_list[[iter_id]]$survey_age_comp$survey1 / rowSums(om_output_list[[iter_id]]$survey_age_comp$survey1) + + for (i in 1:length(c(t(om_cnaa_proportion)))) { + expect_equal(c(t(fims_cnaa_proportion))[i], c(t(om_cnaa_proportion))[i]) + } +}) + +test_that("nll test of fims", { + + result <- setup_and_run_FIMS_without_wrappers( + 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 + ) + + # Set up TMB's computational graph + obj <- result$obj + report <- result$report + + # Calculate standard errors + sdr <- result$sdr + sdr_fixed <- result$sdr_fixed + + # log(R0) + fims_logR0 <- sdr_fixed[36, "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)) + + # recruitment likelihood + # log_devs is of length nyr-1 + rec_nll <- -sum(dnorm( + om_input_list[[iter_id]]$logR.resid[-1], rep(0, om_input_list[[iter_id]]$nyr - 1), + om_input_list[[iter_id]]$logR_sd, TRUE + )) + + # catch and survey index expected likelihoods + 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(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 + )) + index_nll <- index_nll_fleet + index_nll_survey + # age comp likelihoods + fishing_acomp_observed <- em_input_list[[iter_id]]$L.age.obs$fleet1 + fishing_acomp_expected <- om_output_list[[iter_id]]$L.age$fleet1 / rowSums(om_output_list[[iter_id]]$L.age$fleet1) + survey_acomp_observed <- em_input_list[[iter_id]]$survey.age.obs$survey1 + survey_acomp_expected <- om_output_list[[iter_id]]$survey_age_comp$survey1 / rowSums(om_output_list[[iter_id]]$survey_age_comp$survey1) + age_comp_nll_fleet <- age_comp_nll_survey <- 0 + for (y in 1:om_input_list[[iter_id]]$nyr) { + age_comp_nll_fleet <- age_comp_nll_fleet - + dmultinom( + fishing_acomp_observed[y, ] * em_input_list[[iter_id]]$n.L$fleet1, em_input_list[[iter_id]]$n.L$fleet1, + fishing_acomp_expected[y, ], TRUE + ) + + age_comp_nll_survey <- age_comp_nll_survey - + dmultinom( + survey_acomp_observed[y, ] * em_input_list[[iter_id]]$n.survey$survey1, em_input_list[[iter_id]]$n.survey$survey1, + survey_acomp_expected[y, ], TRUE + ) + } + age_comp_nll <- age_comp_nll_fleet + age_comp_nll_survey + + # length comp likelihoods + # TODO: the commented-out code below is not working yet + # fishing_lengthcomp_observed <- em_input_list[[iter_id]]$L.length.obs$fleet1 + # fishing_lengthcomp_expected <- om_output_list[[iter_id]]$L.length$fleet1 / rowSums(om_output_list[[iter_id]]$L.length$fleet1) + # survey_lengthcomp_observed <- em_input_list[[iter_id]]$survey.length.obs$survey1 + # survey_lengthcomp_expected <- om_output_list[[iter_id]]$survey_length_comp$survey1 / rowSums(om_output_list[[iter_id]]$survey_length_comp$survey1) + # lengthcomp_nll_fleet <- lengthcomp_nll_survey <- 0 + # for (y in 1:om_input_list[[iter_id]]$nyr) { + # lengthcomp_nll_fleet <- lengthcomp_nll_fleet - + # dmultinom( + # fishing_lengthcomp_observed[y, ] * em_input_list[[iter_id]]$n.L.lengthcomp$fleet1, em_input_list[[iter_id]]$n.L.lengthcomp$fleet1, + # fishing_lengthcomp_expected[y, ], TRUE + # ) + # + # lengthcomp_nll_survey <- lengthcomp_nll_survey - + # dmultinom( + # survey_lengthcomp_observed[y, ] * em_input_list[[iter_id]]$n.survey.lengthcomp$survey1, em_input_list[[iter_id]]$n.survey.lengthcomp$survey1, + # survey_lengthcomp_expected[y, ], TRUE + # ) + # } + # lengthcomp_nll <- lengthcomp_nll_fleet + lengthcomp_nll_survey + # + # expected_jnll <- rec_nll + index_nll + age_comp_nll + lengthcomp_nll + jnll <- report$jnll + + expect_equal(report$nll_components[1], rec_nll) + expect_equal(report$nll_components[2], index_nll_fleet) + expect_equal(report$nll_components[3], age_comp_nll_fleet) + # expect_equal(report$nll_components[4], lengthcomp_nll_fleet) + expect_equal(report$nll_components[5], index_nll_survey) + expect_equal(report$nll_components[6], age_comp_nll_survey) + # expect_equal(report$nll_components[7], lengthcomp_nll_survey) + # expect_equal(report$jnll, expected_jnll) +}) + +test_that("estimation test of fims", { + + result <- setup_and_run_FIMS_without_wrappers( + iter_id = iter_id, + om_input_list = om_input_list, + om_output_list = om_output_list, + em_input_list = em_input_list, + estimation_mode = TRUE + ) + + # Compare FIMS results with model comparison project OM values + validate_fims( + report = result$report, + sdr = result$sdr, + sdr_report = result$sdr_report, + om_input = om_input_list[[iter_id]], + om_output = om_output_list[[iter_id]], + em_input = em_input_list[[iter_id]] + ) +}) + +test_that("run FIMS with missing values", { + + # Define the NA (missing value) placeholder and the index where it will be inserted + na_value <- -999 + na_index <- 2 + + # Introduce a missing value into the survey observations for the estimation model input + em_input_list[[iter_id]]$surveyB.obs$survey1[na_index] <- na_value + + # Run the FIMS setup and execution function + result <- setup_and_run_FIMS_without_wrappers( + iter_id = iter_id, + om_input_list = om_input_list, + om_output_list = om_output_list, + em_input_list = em_input_list, + estimation_mode = TRUE + ) + + # Validate that the result report is not null + expect_false(is.null(result$report)) + + # Obtain the gradient and Hessian matrix + g <- as.numeric(result$obj$gr(result$opt$par)) + h <- optimHess(result$opt$par, fn = result$obj$fn, gr = result$obj$gr) + result$opt$par <- result$opt$par - solve(h, g) + + # Obtain the maximum absolute gradient to check convergence + # Ensure that the maximum gradient is less than or equal to + # the specified tolerance (0.0001) + max_gradient <- max(abs(result$obj$gr(result$opt$par))) + expect_lte(max_gradient, 0.0001) +}) + +test_that("agecomp in proportion works", { + + # Store the original values of the number of landings observations and + # survey observations + n.L_original <- om_input_list[[iter_id]]$n.L$fleet1 + n.survey_original <- om_input_list[[iter_id]]$n.survey$survey1 + + # Set the number of landings observations and survey observations to 1 + om_input_list[[iter_id]]$n.L$fleet1 <- 1 + om_input_list[[iter_id]]$n.survey$survey1 <- 1 + on.exit(om_input_list[[iter_id]]$n.L$fleet1 <- n.L_original, add = TRUE) + on.exit(om_input_list[[iter_id]]$n.survey$survey1 <- n.survey_original, add = TRUE) + + # Run the FIMS setup and execution function + result <- setup_and_run_FIMS_without_wrappers( + iter_id = iter_id, + om_input_list = om_input_list, + om_output_list = om_output_list, + em_input_list = em_input_list, + estimation_mode = TRUE + ) + + # Compare FIMS results with model comparison project OM values + validate_fims( + report = result$report, + sdr = TMB::sdreport(result$obj), + sdr_report = result$sdr_report, + om_input = om_input_list[[iter_id]], + om_output = om_output_list[[iter_id]], + em_input = em_input_list[[iter_id]] + ) +}) diff --git a/tests/testthat/test-parallel-with-snowfall-with-wrappers.R b/tests/testthat/test-parallel-with-snowfall-with-wrappers.R index f517e57d..198d0ba0 100644 --- a/tests/testthat/test-parallel-with-snowfall-with-wrappers.R +++ b/tests/testthat/test-parallel-with-snowfall-with-wrappers.R @@ -1,76 +1,76 @@ -# # Ensure the latest precompiled version of FIMS is installed in R before -# # running devtools. To do this, either run: -# # - devtools::install() followed by devtools::test(), or -# # - devtools::check() +# Ensure the latest precompiled version of FIMS is installed in R before +# running devtools. To do this, either run: +# - devtools::install() followed by devtools::test(), or +# - devtools::check() -# # Run FIMS in serial and parallel -# # This test demonstrates how to run the FIMS model in both serial and parallel -# # modes. The parallel execution uses {snowfall} to parallelize the tasks across -# # multiple CPU cores. +# Run FIMS in serial and parallel +# This test demonstrates how to run the FIMS model in both serial and parallel +# modes. The parallel execution uses {snowfall} to parallelize the tasks across +# multiple CPU cores. -# # Load the model comparison operating model data from the fixtures folder -# load(test_path("fixtures", "integration_test_data.RData")) +# Load the model comparison operating model data from the fixtures folder +load(test_path("fixtures", "integration_test_data.RData")) -# sim_num <- 10 +sim_num <- 10 -# # Run the FIMS model in serial and record the execution time -# estimation_results_serial <- vector(mode = "list", length = sim_num) +# Run the FIMS model in serial and record the execution time +estimation_results_serial <- vector(mode = "list", length = sim_num) -# for (i in 1:sim_num) { -# estimation_results_serial[[i]] <- setup_and_run_FIMS_with_wrappers( -# iter_id = i, -# om_input_list = om_input_list, -# om_output_list = om_output_list, -# em_input_list = em_input_list, -# estimation_mode = TRUE -# ) -# } +for (i in 1:sim_num) { + estimation_results_serial[[i]] <- setup_and_run_FIMS_with_wrappers( + iter_id = i, + om_input_list = om_input_list, + om_output_list = om_output_list, + em_input_list = em_input_list, + estimation_mode = TRUE + ) +} -# test_that("Run FIMS in parallel using {snowfall}", { -# core_num <- 2 -# snowfall::sfInit(parallel = TRUE, cpus = core_num) +test_that("Run FIMS in parallel using {snowfall}", { + core_num <- 2 + snowfall::sfInit(parallel = TRUE, cpus = core_num) -# snowfall::sfLibrary(FIMS) -# results_parallel <- snowfall::sfLapply( -# 1:sim_num, -# setup_and_run_FIMS_with_wrappers, -# om_input_list, -# om_output_list, -# em_input_list, -# TRUE -# ) + snowfall::sfLibrary(FIMS) + results_parallel <- snowfall::sfLapply( + 1:sim_num, + setup_and_run_FIMS_with_wrappers, + om_input_list, + om_output_list, + em_input_list, + TRUE + ) -# snowfall::sfStop() + snowfall::sfStop() -# # Comparison of results: -# # Verify that SSB values from both runs are equivalent. -# expect_setequal( -# purrr::map( -# results_parallel, -# \(x) x@estimates[x@estimates$name == "SSB", "value"] -# ), -# purrr::map( -# estimation_results_serial, -# \(x) x@estimates[x@estimates$name == "SSB", "value"] -# ) -# ) + # Comparison of results: + # Verify that SSB values from both runs are equivalent. + expect_setequal( + purrr::map( + results_parallel, + \(x) x@estimates[x@estimates$name == "SSB", "value"] + ), + purrr::map( + estimation_results_serial, + \(x) x@estimates[x@estimates$name == "SSB", "value"] + ) + ) -# # Verify that parameter values from both runs are equivalent. -# expect_setequal( -# purrr::map( -# results_parallel, -# \(x) x@estimates[x@estimates$name == "p", "value"] -# ), -# purrr::map( -# estimation_results_serial, -# \(x) x@estimates[x@estimates$name == "p", "value"] -# ) -# ) + # Verify that parameter values from both runs are equivalent. + expect_setequal( + purrr::map( + results_parallel, + \(x) x@estimates[x@estimates$name == "p", "value"] + ), + purrr::map( + estimation_results_serial, + \(x) x@estimates[x@estimates$name == "p", "value"] + ) + ) -# # Verify that total NLL values from both runs are equivalent. -# expect_equal( -# purrr::map(results_parallel, \(x) x@report[["jnll"]]), -# purrr::map(estimation_results_serial, \(x) x@report[["jnll"]]) -# ) + # Verify that total NLL values from both runs are equivalent. + expect_equal( + purrr::map(results_parallel, \(x) x@report[["jnll"]]), + purrr::map(estimation_results_serial, \(x) x@report[["jnll"]]) + ) -# }) +}) diff --git a/vignettes/fims-demo.Rmd b/vignettes/fims-demo.Rmd index f47c874b..5311be19 100644 --- a/vignettes/fims-demo.Rmd +++ b/vignettes/fims-demo.Rmd @@ -70,23 +70,33 @@ fleet1 <- survey1 <- list( selectivity = list(form = "LogisticSelectivity"), data_distribution = c( Index = "DlnormDistribution", - AgeComp = "DmultinomDistribution" + AgeComp = "DmultinomDistribution", + LengthComp = "DmultinomDistribution" ) ) # Create default parameters default_parameters <- fims_frame |> create_default_parameters( - fleets = list(fleet1 = fleet1, survey1 = survey1), - recruitment = list( - form = "BevertonHoltRecruitment", - process_distribution = c(log_devs = "DnormDistribution") - ), - growth = list(form = "EWAAgrowth"), - maturity = list(form = "LogisticMaturity") + fleets = list(fleet1 = fleet1, survey1 = survey1) ) ``` -We may update default parameters to experiment with different parameter values. Here, we adjust fishing mortality, selectivity, maturity, and population parameters using `update_parameters()`. + +create_default_parameters() includes default settings for recruitment, growth, and maturity modules. For example: + + - "BevertonHoltRecruitment" for the recruitment module + - "DnormDistribution" for recruitment deviations (log_devs) + - "EWAAgrowth" for the growth module + - "LogisticMaturity" for maturity module + +The argument names and their corresponding default values of the create_default_parameters() function can be displayed using `args(FIMS::create_default_parameters)`. +```{r display-create-default-parameters-default-values} +args(FIMS::create_default_parameters) +``` + +We may update the default parameters to experiment with different values. Before modifying the parameter values, we can print the `default_parameters`, then modify a specific element or a list of elements. The dimensions of the updated elements should match the original dimensions from the `default_parameters`. + +Here, we adjust fishing mortality, selectivity, maturity, and population parameters using `update_parameters()`. ```{r modify-default-parameters} parameters <- default_parameters |> update_parameters( @@ -200,6 +210,51 @@ low_slope_fit <- parameters_low_slope |> initialize_fims(data = fims_frame) |> fit_fims(optimize = TRUE) +clear() + +# Fitting models with only age or only length data +# Fitting With only age data +# Define fleet and survey with age-specific data distribution +fleet1 <- survey1 <- list( + selectivity = list(form = "LogisticSelectivity"), + data_distribution = c( + Index = "DlnormDistribution", + AgeComp = "DmultinomDistribution" + ) +) + +# Create default parameters, update with modified values, initialize FIMS, +# and fit the model +age_only_fit <- fims_frame |> + create_default_parameters( + fleets = list(fleet1 = fleet1, survey1 = survey1) + ) |> + update_parameters(modified_parameters = parameters$parameters) |> + initialize_fims(data = fims_frame) |> + fit_fims(optimize = TRUE) + +clear() + +# Fitting with only length data +# Define fleet and survey with length-specific data distribution +fleet1 <- survey1 <- list( + selectivity = list(form = "LogisticSelectivity"), + data_distribution = c( + Index = "DlnormDistribution", + LengthComp = "DmultinomDistribution" + ) +) + +# Create default parameters, update with modified values, initialize FIMS, +# and fit the model +length_only_fit <- fims_frame |> + create_default_parameters( + fleets = list(fleet1 = fleet1, survey1 = survey1) + ) |> + update_parameters(modified_parameters = parameters$parameters) |> + initialize_fims(data = fims_frame) |> + fit_fims(optimize = TRUE) + clear() ```