diff --git a/tests/testthat/helper-integration-tests-setup.R b/tests/testthat/helper-integration-tests-setup.R index b64568b2..4782f8b8 100644 --- a/tests/testthat/helper-integration-tests-setup.R +++ b/tests/testthat/helper-integration-tests-setup.R @@ -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. diff --git a/tests/testthat/test-integration-fims-estimation-with-wrappers.R b/tests/testthat/test-integration-fims-estimation-with-wrappers.R index 95b26d5e..81497bcf 100644 --- a/tests/testthat/test-integration-fims-estimation-with-wrappers.R +++ b/tests/testthat/test-integration-fims-estimation-with-wrappers.R @@ -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", { diff --git a/tests/testthat/test-integration-fims-estimation-without-wrappers.R b/tests/testthat/test-integration-fims-estimation-without-wrappers.R index 6c2d243e..d8c1f457 100644 --- a/tests/testthat/test-integration-fims-estimation-without-wrappers.R +++ b/tests/testthat/test-integration-fims-estimation-without-wrappers.R @@ -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) 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", {