Skip to content

Commit

Permalink
make FIMS run without age composition data
Browse files Browse the repository at this point in the history
* add if statement in functions from R/initialize_modules.R to create a FIMS model without age comp data
* add an R test to run a FIMS model without age comp data
  • Loading branch information
Bai-Li-NOAA committed Dec 20, 2024
1 parent af49779 commit 3326ddc
Show file tree
Hide file tree
Showing 2 changed files with 120 additions and 19 deletions.
45 changes: 27 additions & 18 deletions R/initialize_modules.R
Original file line number Diff line number Diff line change
Expand Up @@ -408,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,
Expand All @@ -423,12 +417,17 @@ 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()

if ("age" %in% fleet_types &
"AgeComp" %in% names(parameters[["modules"]][["fleets"]][[fleet_name]][["data_distribution"]])) {
module$SetObservedAgeCompData(linked_ids["age_comp"])
}

if ("length" %in% fleet_types &
"LengthComp" %in% names(parameters[["modules"]][["fleets"]][[fleet_name]][["data_distribution"]])) {
module$SetObservedLengthCompData(linked_ids["length_comp"])
Expand Down Expand Up @@ -616,14 +615,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()
)

Expand All @@ -632,6 +625,19 @@ initialize_fims <- function(parameters, data) {
dplyr::pull(type) |>
unique()

if ("age" %in% fleet_types &
"AgeComp" %in% names(parameters[["modules"]][["fleets"]][[fleet_names[i]]][["data_distribution"]])) {
fleet_age_comp[[i]] <- initialize_age_comp(
data = data,
fleet_name = fleet_names[i]
)

fleet_module_ids <- c(
fleet_module_ids,
c(age_comp = fleet_age_comp[[i]]$get_id())
)
}

if ("length" %in% fleet_types &
"LengthComp" %in% names(parameters[["modules"]][["fleets"]][[fleet_names[i]]][["data_distribution"]])) {
fleet_length_comp[[i]] <- initialize_length_comp(
Expand Down Expand Up @@ -683,11 +689,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"]])) {
Expand Down
94 changes: 93 additions & 1 deletion tests/testthat/test-integration-fims-estimation-with-wrappers.R
Original file line number Diff line number Diff line change
Expand Up @@ -231,7 +231,7 @@ test_that("estimation test of fims using wrapper functions", {
)
})

test_that("estimation test with length comp using 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]]
Expand Down Expand Up @@ -322,3 +322,95 @@ test_that("estimation test with length comp using wrappers",{
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 = 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
)
})

0 comments on commit 3326ddc

Please sign in to comment.