Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix length tests #739

Merged
merged 1 commit into from
Jan 29, 2025
Merged

Fix length tests #739

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 18 additions & 0 deletions tests/testthat/helper-integration-tests-setup.R
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,3 +1,21 @@
#' FIMS dmultinom()
#' This function matches the dmultinom() function in TMB and differs from R
#' by NOT rounding obs to the nearest integer. The function is evaluated in
#' log space and returns the log probability mass function.
#'
#' @param x A vector of length K of numeric values.
#'
#' @param p A numeric non-negative vector of length K, specifying the probability
#' for the K classes; must sum 1.
#'
#' @return The log of the probability mass function for the multinomal.
FIMS_dmultinom <- function(x, p){
xp1 <- x + 1
log_pmf = lgamma(sum(x) + 1) - sum(lgamma(xp1)) + sum(x*log(p))
return(log_pmf)
}


#' Set Up and Run FIMS Model without using wrapper functions
#'
#' This function sets up and runs the FIMS for a given iteration.
Expand Down
54 changes: 27 additions & 27 deletions tests/testthat/test-integration-fims-estimation-with-wrappers.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -199,39 +199,39 @@ test_that("nll test of fims", {
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"]]
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"]]) {
# test using FIMS_dmultinom which matches the TMB dmultinom calculation and differs from R
# by NOT rounding obs to the nearest integer.
lengthcomp_nll_fleet <- lengthcomp_nll_fleet -
FIMS_dmultinom(
fishing_lengthcomp_observed[y, ] * em_input_list[[iter_id]][["n.L.lengthcomp"]][["fleet1"]],
fishing_lengthcomp_expected[y, ]
)

lengthcomp_nll_survey <- lengthcomp_nll_survey -
FIMS_dmultinom(
survey_lengthcomp_observed[y, ] * em_input_list[[iter_id]][["n.survey.lengthcomp"]][["survey1"]],
survey_lengthcomp_expected[y, ]
)
}
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"]][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)
expect_equal(report[["nll_components"]][7], lengthcomp_nll_survey)
expect_equal(jnll, expected_jnll)
})

test_that("estimation test of fims using wrapper functions", {
Expand Down
51 changes: 26 additions & 25 deletions tests/testthat/test-integration-fims-estimation-without-wrappers.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -210,38 +210,39 @@ test_that("nll test of fims", {
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
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"]]) {
# test using FIMS_dmultinom which matches the TMB dmultinom calculation and differs from R
# by NOT rounding obs to the nearest integer.
lengthcomp_nll_fleet <- lengthcomp_nll_fleet -
FIMS_dmultinom(
fishing_lengthcomp_observed[y, ] * em_input_list[[iter_id]][["n.L.lengthcomp"]][["fleet1"]],
fishing_lengthcomp_expected[y, ]
)

lengthcomp_nll_survey <- lengthcomp_nll_survey -
FIMS_dmultinom(
survey_lengthcomp_observed[y, ] * em_input_list[[iter_id]][["n.survey.lengthcomp"]][["survey1"]],
survey_lengthcomp_expected[y, ]
)
}
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)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Non-blocking suggestion for this PR: Could we assign descriptive names to individual nll_components? That way, we won't have to guess that the first element corresponds to recruitment nll and the second to fleet1 index nll.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Great idea @Bai-Li-NOAA. I am just not entirely sure how to get this done because the elements present will change given the model structure. I am sure it is possible but we might have to think hard about how to do it correctly.

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"]][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)
expect_equal(report[["nll_components"]][7], lengthcomp_nll_survey)
expect_equal(report[["jnll"]], expected_jnll)
})

test_that("estimation test of fims", {
Expand Down
Loading