From a98c48e140b9ea15db59dea9f6a52da12201e047 Mon Sep 17 00:00:00 2001
From: Andrea-Havron-NOAA <andrea.havron@noaa.gov>
Date: Thu, 8 Aug 2024 10:31:32 -0700
Subject: [PATCH] fix tests

---
 .../test-integration-fims-estimation.R        | 1289 +++++++++--------
 tests/testthat/test-rcpp-tmb-distributions.R  |  103 +-
 2 files changed, 676 insertions(+), 716 deletions(-)

diff --git a/tests/testthat/test-integration-fims-estimation.R b/tests/testthat/test-integration-fims-estimation.R
index 6f26cf62b..91d97cc88 100644
--- a/tests/testthat/test-integration-fims-estimation.R
+++ b/tests/testthat/test-integration-fims-estimation.R
@@ -56,13 +56,13 @@ setup_fims <- function(om_input, om_output, em_input) {
   # recruit deviations should enter the model in normal space.
   # The log is taken in the likelihood calculations
   # alternative setting: recruitment$log_devs <- rep(0, length(om_input$logR.resid))
-  test_env$recruitment$log_devs <- methods::new(test_env$ParameterVector, om_input$logR.resid[-1], length(om_input$logR.resid[-1]))
+  test_env$recruitment$log_devs <- methods::new(test_env$fims$ParameterVector, om_input$logR.resid[-1], length(om_input$logR.resid[-1]))
 
-  test_env$recruitment_distribution <- new(test_env$TMBDnormDistribution)
+  test_env$recruitment_distribution <- new(test_env$fims$TMBDnormDistribution)
   # set up logR_sd using the normal log_sd parameter
   # logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is
   # taken before the likelihood calculation
-  test_env$recruitment_distribution$log_sd <- new(test_env$ParameterVector, 1)
+  test_env$recruitment_distribution$log_sd <- new(test_env$fims$ParameterVector, 1)
   test_env$recruitment_distribution$log_sd[1]$value <- log(om_input$logR_sd)
   test_env$recruitment_distribution$log_sd[1]$estimated = FALSE
   test_env$recruitment_distribution$set_distribution_links("random_effects", test_env$recruitment$log_devs$get_id())
@@ -113,7 +113,7 @@ setup_fims <- function(om_input, om_output, em_input) {
   test_env$fishing_fleet <- new(test_env$fims$Fleet)
   test_env$fishing_fleet$nages <- om_input$nages
   test_env$fishing_fleet$nyears <- om_input$nyr
-  test_env$fishing_fleet$log_Fmort <- new(test_env$ParameterVector, log(om_output$f), om_input$nyr)
+  test_env$fishing_fleet$log_Fmort <- new(test_env$fims$ParameterVector, log(om_output$f), om_input$nyr)
   test_env$fishing_fleet$log_Fmort$set_all_estimable(TRUE)
   #test_env$fishing_fleet$random_F <- FALSE
   test_env$fishing_fleet$log_q <- log(1.0)
@@ -121,20 +121,20 @@ setup_fims <- function(om_input, om_output, em_input) {
   test_env$fishing_fleet$random_q <- FALSE
   test_env$fishing_fleet$SetSelectivity(test_env$fishing_fleet_selectivity$get_id())
 
-  # Set up fishery index data using the lognormal
-  test_env$fishing_fleet_index_distribution <- methods::new(test_env$TMBDlnormDistribution)
-  #lognormal observation error transformed on the log scale
-  test_env$fishing_fleet_index_distribution$log_logsd <- new(test_env$ParameterVector, om_input$nyr)
+  # Set up fishery index data using the normal
+  test_env$fishing_fleet_index_distribution <- methods::new(test_env$fims$TMBDnormDistribution)
+  #normal observation error transformed on the log scale
+  test_env$fishing_fleet_index_distribution$log_sd <- new(test_env$fims$ParameterVector, om_input$nyr)
   for(y in 1:om_input$nyr){
-    test_env$fishing_fleet_index_distribution$log_logsd[y]$value <- log(sqrt(log(em_input$cv.L$fleet1^2 + 1)))
+    test_env$fishing_fleet_index_distribution$log_sd[y]$value <- log(sqrt(log(em_input$cv.L$fleet1^2 + 1)))
   }
-  test_env$fishing_fleet_index_distribution$log_logsd$set_all_estimable(FALSE)
+  test_env$fishing_fleet_index_distribution$log_sd$set_all_estimable(FALSE)
   # Set Data using the IDs from the modules defined above
   test_env$fishing_fleet_index_distribution$set_observed_data(test_env$fishing_fleet_index$get_id())
   test_env$fishing_fleet_index_distribution$set_distribution_links("data", test_env$fishing_fleet$log_expected_index$get_id())
 
   # Set up fishery age composition data using the multinomial
-  test_env$fishing_fleet_agecomp_distribution <- methods::new(test_env$TMBDmultinomDistribution)
+  test_env$fishing_fleet_agecomp_distribution <- methods::new(test_env$fims$TMBDmultinomDistribution)
   test_env$fishing_fleet_agecomp_distribution$set_observed_data(test_env$fishing_fleet_age_comp$get_id())
   test_env$fishing_fleet_agecomp_distribution$set_distribution_links("data", test_env$fishing_fleet$proportion_catch_numbers_at_age$get_id())
 
@@ -161,33 +161,33 @@ setup_fims <- function(om_input, om_output, em_input) {
   test_env$survey_fleet$random_q <- FALSE
   test_env$survey_fleet$SetSelectivity(test_env$survey_fleet_selectivity$get_id())
 
-  # Set up survey index data using the lognormal
-  test_env$survey_fleet_index_distribution <- methods::new(test_env$TMBDlnormDistribution)
-  #lognormal observation error transformed on the log scale
+  # Set up survey index data using the normal
+  test_env$survey_fleet_index_distribution <- methods::new(test_env$fims$TMBDnormDistribution)
+  #normal observation error transformed on the log scale
   # sd = sqrt(log(cv^2 + 1)), sd is log transformed
-  test_env$survey_fleet_index_distribution$log_logsd <- new(test_env$ParameterVector, om_input$nyr)
+  test_env$survey_fleet_index_distribution$log_sd <- new(test_env$fims$ParameterVector, om_input$nyr)
   for(y in 1:om_input$nyr){
-    test_env$survey_fleet_index_distribution$log_logsd[y]$value <- log(sqrt(log(em_input$cv.survey$survey1^2 + 1)))
+    test_env$survey_fleet_index_distribution$log_sd[y]$value <- log(sqrt(log(em_input$cv.survey$survey1^2 + 1)))
   }
-  test_env$survey_fleet_index_distribution$log_logsd$set_all_estimable(FALSE)
+  test_env$survey_fleet_index_distribution$log_sd$set_all_estimable(FALSE)
   # Set Data using the IDs from the modules defined above
   test_env$survey_fleet_index_distribution$set_observed_data(test_env$survey_fleet_index$get_id())
   test_env$survey_fleet_index_distribution$set_distribution_links("data", test_env$survey_fleet$log_expected_index$get_id())
 
   # Age composition data
 
-  test_env$survey_fleet_agecomp_distribution <- methods::new(test_env$TMBDmultinomDistribution)
+  test_env$survey_fleet_agecomp_distribution <- methods::new(test_env$fims$TMBDmultinomDistribution)
   test_env$survey_fleet_agecomp_distribution$set_observed_data(test_env$survey_fleet_age_comp$get_id())
   test_env$survey_fleet_agecomp_distribution$set_distribution_links("data", test_env$survey_fleet$proportion_catch_numbers_at_age$get_id())
 
 
   # Population
   test_env$population <- new(test_env$fims$Population)
-  test_env$population$log_M <- methods::new(test_env$ParameterVector,
+  test_env$population$log_M <- methods::new(test_env$fims$ParameterVector,
     rep(log(om_input$M.age[1]), om_input$nyr * om_input$nages),
      om_input$nyr * om_input$nages)
   test_env$population$log_M$set_all_estimable(FALSE)
-  test_env$population$log_init_naa <- methods::new(test_env$ParameterVector,
+  test_env$population$log_init_naa <- methods::new(test_env$fims$ParameterVector,
     log(om_output$N.age[1, ]), om_input$nages)
   test_env$population$log_init_naa$set_all_estimable(TRUE)
   test_env$population$nages <- om_input$nages
@@ -202,442 +202,447 @@ setup_fims <- function(om_input, om_output, em_input) {
   # end of setup_fims function, returning test_env
   return(test_env)
 }
-#
-# test_that("deterministic test of fims", {
-#   # run function defined above to set up the test environment
-#   deterministic_env <- setup_fims(
-#     om_input = om_input,
-#     om_output = om_output,
-#     em_input = em_input
-#   )
-#   # Set-up
-#   deterministic_env$fims$CreateTMBModel()
-#   # CreateTMBModel calls a function in information that loops
-#   # over all the populations and fleets and sets all the pointers
-#   parameters <- list(p = deterministic_env$fims$get_fixed())
-#   # get_fixed function is an Rcpp function that loops over all Rcpp
-#   # modules and returned a vector of parameters being estimated
-#
-#   # Set up TMB's computational graph
-#   obj <- MakeADFun(data = list(), parameters, DLL = "FIMS")
-#
-#   # Calculate standard errors
-#   sdr <- TMB::sdreport(obj)
-#   sdr_fixed <- summary(sdr, "fixed")
-#
-#   # Call report using deterministic parameter values
-#   # obj$report() requires parameter list to avoid errors
-#   report <- obj$report(obj$par)
-#
-#   # Compare log(R0) to true value
-#   fims_logR0 <- sdr_fixed[1, "Estimate"]
-#   expect_gt(fims_logR0, 0.0)
-#   expect_equal(fims_logR0, log(om_input$R0))
-#
-#   # Compare numbers at age to true value
-#   for (i in 1:length(c(t(om_output$N.age)))) {
-#     expect_equal(report$naa[[1]][i], c(t(om_output$N.age))[i])
-#   }
-#
-#   # Compare biomass to true value
-#   for (i in 1:length(om_output$biomass.mt)) {
-#     expect_equal(report$biomass[[1]][i], om_output$biomass.mt[i])
-#   }
-#
-#   # Compare spawning biomass to true value
-#   for (i in 1:length(om_output$SSB)) {
-#     expect_equal(report$ssb[[1]][i], om_output$SSB[i])
-#   }
-#
-#   # Compare recruitment to true value
-#   fims_naa <- matrix(report$naa[[1]][1:(om_input$nyr * om_input$nages)],
-#     nrow = om_input$nyr, byrow = TRUE
-#   )
-#
-#   # loop over years to compare recruitment by year
-#   for (i in 1:om_input$nyr) {
-#     expect_equal(fims_naa[i, 1], om_output$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$nyr, 1],
-#     report$recruitment[[1]][1:om_input$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$nyr) {
-#     expect_equal(report$recruitment[[1]][i], om_output$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$logR.resid[-1])
-#
-#   # F (fixed at initial "true" values)
-#   expect_equal(report$F_mort[[1]], om_output$f)
-#
-#   # Expected catch
-#   fims_index <- report$exp_index
-#   for (i in 1:length(om_output$L.mt$fleet1)) {
-#     expect_equal(fims_index[[1]][i], om_output$L.mt$fleet1[i])
-#   }
-#
-#   # Expect small relative error for deterministic test
-#   fims_object_are <- rep(0, length(em_input$L.obs$fleet1))
-#   for (i in 1:length(em_input$L.obs$fleet1)) {
-#     fims_object_are[i] <- abs(fims_index[[1]][i] - em_input$L.obs$fleet1[i]) / em_input$L.obs$fleet1[i]
-#   }
-#
-#   # Expect 95% of relative error to be within 2*cv
-#   expect_lte(sum(fims_object_are > om_input$cv.L$fleet1 * 2.0), length(em_input$L.obs$fleet1) * 0.05)
-#
-#   # Compare expected catch number at age to true values
-#   for (i in 1:length(c(t(om_output$L.age$fleet1)))) {
-#     expect_equal(report$cnaa[[1]][i], c(t(om_output$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$nyr * om_input$nages)],
-#     nrow = om_input$nyr, byrow = TRUE
-#   )
-#   fims_cnaa_proportion <- fims_cnaa / rowSums(fims_cnaa)
-#   om_cnaa_proportion <- om_output$L.age$fleet1 / rowSums(om_output$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$nyr * om_input$nages)], nrow = om_input$nyr, byrow = TRUE)
-#   expect_equal(fims_index[[2]], apply(cwaa, 1, sum) * om_output$survey_q$survey1)
-#
-#   for (i in 1:length(om_output$survey_index_biomass$survey1)) {
-#     expect_equal(fims_index[[2]][i], om_output$survey_index_biomass$survey1[i])
-#   }
-#
-#   fims_object_are <- rep(0, length(em_input$surveyB.obs$survey1))
-#   for (i in 1:length(em_input$survey.obs$survey1)) {
-#     fims_object_are[i] <- abs(fims_index[[2]][i] - em_input$surveyB.obs$survey1[i]) / em_input$surveyB.obs$survey1[i]
-#   }
-#   # Expect 95% of relative error to be within 2*cv
-#   expect_lte(sum(fims_object_are > om_input$cv.survey$survey1 * 2.0), length(em_input$surveyB.obs$survey1) * 0.05)
-#
-#   # Expected catch number at age in proportion
-#   fims_cnaa <- matrix(report$cnaa[[2]][1:(om_input$nyr * om_input$nages)],
-#     nrow = om_input$nyr, byrow = TRUE
-#   )
-#
-#   for (i in 1:length(c(t(om_output$survey_age_comp$survey1)))) {
-#     expect_equal(report$cnaa[[2]][i], c(t(om_output$survey_age_comp$survey1))[i])
-#   }
-#
-#   fims_cnaa_proportion <- fims_cnaa / rowSums(fims_cnaa)
-#   om_cnaa_proportion <- om_output$survey_age_comp$survey1 / rowSums(om_output$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])
-#   }
-#   # clear memory
-#   deterministic_env$fims$clear()
-# })
-#
-# # test likelihood components once new logging system is working
-# # test_that("nll test of fims", {
-# #   nll_env <- setup_fims(
-# #     om_input = om_input,
-# #     om_output = om_output,
-# #     em_input = em_input
-# #   )
-# #   # Set-up TMB
-# #   nll_env$fims$CreateTMBModel()
-# #   parameters <- list(p = nll_env$fims$get_fixed())
-# #   par_list <- 1:length(parameters[[1]])
-# #   par_list[2:length(par_list)] <- NA
-# #   map <- list(p = factor(par_list))
-# #
-# #   obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS", map = map)
-# #
-# #   sdr <- TMB::sdreport(obj)
-# #   sdr_fixed <- summary(sdr, "fixed")
-# #
-# #   # log(R0)
-# #   fims_logR0 <- sdr_fixed[1, "Estimate"]
-# #   # expect_lte(abs(fims_logR0 - log(om_input$R0)) / log(om_input$R0), 0.0001)
-# #   expect_equal(fims_logR0, log(om_input$R0))
-# #
-# #   # Call report using deterministic parameter values
-# #   # obj$report() requires parameter list to avoid errors
-# #   report <- obj$report(obj$par)
-# #   obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS", map = map)
-# #   jnll <- obj$fn()
-# #
-# #
-# #   # recruitment likelihood
-# #   # log_devs is of length nyr-1
-# #   rec_nll <- -sum(dnorm(
-# #     nll_env$recruitment$log_devs, rep(0, om_input$nyr - 1),
-# #     om_input$logR_sd, TRUE
-# #   ))
-# #
-# #   # catch and survey index expected likelihoods
-# #   index_nll_fleet <- -sum(dnorm(
-# #     log(nll_env$catch),
-# #     log(om_output$L.mt$fleet1),
-# #     sqrt(log(em_input$cv.L$fleet1^2 + 1)), TRUE
-# #   ))
-# #   index_nll_survey <- -sum(dnorm(
-# #     log(nll_env$survey_index),
-# #     log(om_output$survey_index_biomass$survey1),
-# #     sqrt(log(em_input$cv.survey$survey1^2 + 1)), TRUE
-# #   ))
-# #   index_nll <- index_nll_fleet + index_nll_survey
-# #   # age comp likelihoods
-# #   fishing_acomp_observed <- em_input$L.age.obs$fleet1
-# #   fishing_acomp_expected <- om_output$L.age$fleet1 / rowSums(om_output$L.age$fleet1)
-# #   survey_acomp_observed <- em_input$survey.age.obs$survey1
-# #   survey_acomp_expected <- om_output$survey_age_comp$survey1 / rowSums(om_output$survey_age_comp$survey1)
-# #   age_comp_nll_fleet <- age_comp_nll_survey <- 0
-# #   for (y in 1:om_input$nyr) {
-# #     age_comp_nll_fleet <- age_comp_nll_fleet -
-# #       dmultinom(
-# #         fishing_acomp_observed[y, ] * em_input$n.L$fleet1, em_input$n.L$fleet1,
-# #         fishing_acomp_expected[y, ], TRUE
-# #       )
-# #
-# #     age_comp_nll_survey <- age_comp_nll_survey -
-# #       dmultinom(
-# #         survey_acomp_observed[y, ] * em_input$n.survey$survey1, em_input$n.survey$survey1,
-# #         survey_acomp_expected[y, ], TRUE
-# #       )
-# #   }
-# #   age_comp_nll <- age_comp_nll_fleet + age_comp_nll_survey
-# #   expected_jnll <- rec_nll + index_nll + age_comp_nll
-# #
-# #   expect_equal(report$rec_nll, rec_nll)
-# #   expect_equal(report$age_comp_nll, age_comp_nll)
-# #   expect_equal(report$index_nll, index_nll)
-# #   expect_equal(jnll, expected_jnll)
-# #
-# #   nll_env$fims$clear()
-# # })
-#
-# test_that("estimation test of fims", {
-#   estimation_env <- setup_fims(
-#     om_input = om_input,
-#     om_output = om_output,
-#     em_input = em_input
-#   )
-#
-#   # Set-up TMB
-#   estimation_env$fims$CreateTMBModel()
-#   # Create parameter list from Rcpp modules
-#   parameters <- list(p = estimation_env$fims$get_fixed())
-#   obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS")
-#
-#   opt <- with(obj, optim(par, fn, gr,
-#     method = "BFGS",
-#     control = list(maxit = 1000000, reltol = 1e-15)
-#   ))
-#
-#   # Call report using MLE parameter values
-#   # obj$report() requires parameter list to avoid errors
-#   report <- obj$report(obj$env$last.par.best)
-#   sdr <- TMB::sdreport(obj)
-#   sdr_report <- summary(sdr, "report")
-#   # Numbers at age
-#   # Estimates and SE for NAA
-#   sdr_naa <- sdr_report[which(rownames(sdr_report) == "NAA"), ]
-#   naa_are <- rep(0, length(c(t(om_output$N.age))))
-#   for (i in 1:length(c(t(om_output$N.age)))) {
-#     naa_are[i] <- abs(sdr_naa[i, 1] - c(t(om_output$N.age))[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of NAA
-#   expect_lte(
-#     sum(naa_are > qnorm(.975) * sdr_naa[1:length(c(t(om_output$N.age))), 2]),
-#     0.05 * length(c(t(om_output$N.age)))
-#   )
-#
-#   # Biomass
-#   sdr_biomass <- sdr_report[which(rownames(sdr_report) == "Biomass"), ]
-#   biomass_are <- rep(0, length(om_output$biomass.mt))
-#   for (i in 1:length(om_output$biomass.mt)) {
-#     biomass_are[i] <- abs(sdr_biomass[i, 1] - om_output$biomass.mt[i]) # / om_output$biomass.mt[i]
-#     # expect_lte(biomass_are[i], 0.15)
-#   }
-#   expect_lte(
-#     sum(biomass_are > qnorm(.975) * sdr_biomass[1:length(om_output$biomass.mt), 2]),
-#     0.05 * length(om_output$biomass.mt)
-#   )
-#
-#   # Spawning biomass
-#   sdr_sb <- sdr_report[which(rownames(sdr_report) == "SSB"), ]
-#   sb_are <- rep(0, length(om_output$SSB))
-#   for (i in 1:length(om_output$SSB)) {
-#     sb_are[i] <- abs(sdr_sb[i, 1] - om_output$SSB[i]) # / om_output$SSB[i]
-#     # expect_lte(sb_are[i], 0.15)
-#   }
-#   expect_lte(
-#     sum(sb_are > qnorm(.975) * sdr_sb[1:length(om_output$SSB), 2]),
-#     0.05 * length(om_output$SSB)
-#   )
-#
-#   # Recruitment
-#   fims_naa <- matrix(report$naa[[1]][1:(om_input$nyr * om_input$nages)],
-#     nrow = om_input$nyr, byrow = TRUE
-#   )
-#   sdr_naa1_vec <- sdr_report[which(rownames(sdr_report) == "NAA"), 2]
-#   sdr_naa1 <- sdr_naa1_vec[seq(1, om_input$nyr * om_input$nages, by = om_input$nages)]
-#   fims_naa1_are <- rep(0, om_input$nyr)
-#   for (i in 1:om_input$nyr) {
-#     fims_naa1_are[i] <- abs(fims_naa[i, 1] - om_output$N.age[i, 1]) # /
-#     # om_output$N.age[i, 1]
-#     # expect_lte(fims_naa1_are[i], 0.25)
-#   }
-#   expect_lte(
-#     sum(fims_naa1_are > qnorm(.975) * sdr_naa1[1:length(om_output$SSB)]),
-#     0.05 * length(om_output$SSB)
-#   )
-#
-#   expect_equal(
-#     fims_naa[, 1],
-#     report$recruitment[[1]][1:om_input$nyr]
-#   )
-#
-#   # recruitment log deviations
-#   # the initial value of om_input$logR.resid is dropped from the model
-#   sdr_rdev <- sdr_report[which(rownames(sdr_report) == "LogRecDev"), ]
-#   rdev_are <- rep(0, length(om_input$logR.resid) - 1)
-#
-#   for (i in 1:(length(report$log_recruit_dev[[1]]) - 1)) {
-#     rdev_are[i] <- abs(report$log_recruit_dev[[1]][i] - om_input$logR.resid[i + 1]) # /
-#     #   exp(om_input$logR.resid[i])
-#     # expect_lte(rdev_are[i], 1) # 1
-#   }
-#   expect_lte(
-#     sum(rdev_are > qnorm(.975) * sdr_rdev[1:length(om_input$logR.resid) - 1, 2]),
-#     0.05 * length(om_input$logR.resid)
-#   )
-#
-#   # F (needs to be updated when std.error is available)
-#   sdr_F <- sdr_report[which(rownames(sdr_report) == "FMort"), ]
-#   f_are <- rep(0, length(om_output$f))
-#   for (i in 1:length(om_output$f)) {
-#     f_are[i] <- abs(sdr_F[i, 1] - om_output$f[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of Fmort
-#   expect_lte(
-#     sum(f_are > qnorm(.975) * sdr_F[1:length(om_output$f), 2]),
-#     0.05 * length(om_output$f)
-#   )
-#
-#   # Expected fishery catch and survey index
-#   fims_index <- sdr_report[which(rownames(sdr_report) == "ExpectedIndex"), ]
-#   fims_catch <- fims_index[1:om_input$nyr, ]
-#   fims_survey <- fims_index[(om_input$nyr + 1):(om_input$nyr * 2), ]
-#
-#   # Expected fishery catch - om_output
-#   catch_are <- rep(0, length(om_output$L.mt$fleet1))
-#   for (i in 1:length(om_output$L.mt$fleet1)) {
-#     catch_are[i] <- abs(fims_catch[i, 1] - om_output$L.mt$fleet1[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of fishery catch
-#   expect_lte(
-#     sum(catch_are > qnorm(.975) * fims_catch[, 2]),
-#     0.05 * length(om_output$L.mt$fleet1)
-#   )
-#
-#   # Expected fishery catch - em_input
-#   catch_are <- rep(0, length(em_input$L.obs$fleet1))
-#   for (i in 1:length(em_input$L.obs$fleet1)) {
-#     catch_are[i] <- abs(fims_catch[i, 1] - em_input$L.obs$fleet1[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of fishery catch
-#   expect_lte(
-#     sum(catch_are > qnorm(.975) * fims_catch[, 2]),
-#     0.05 * length(em_input$L.obs$fleet1)
-#   )
-#
-#
-#   # Expected fishery catch number at age
-#   sdr_cnaa <- sdr_report[which(rownames(sdr_report) == "CNAA"), ]
-#   cnaa_are <- rep(0, length(c(t(om_output$L.age$fleet1))))
-#   for (i in 1:length(c(t(om_output$L.age$fleet1)))) {
-#     cnaa_are[i] <- abs(sdr_cnaa[i, 1] - c(t(om_output$L.age$fleet1))[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of CNAA
-#   expect_lte(
-#     sum(cnaa_are > qnorm(.975) * sdr_cnaa[, 2]),
-#     0.05 * length(c(t(om_output$L.age$fleet1)))
-#   )
-#
-#
-#   # Expected catch number at age in proportion
-#   # fims_cnaa <- matrix(report$cnaa[1:(om_input$nyr*om_input$nages), 1],
-#   #                     nrow = om_input$nyr, byrow = TRUE)
-#   # fims_cnaa_proportion <- fims_cnaa/rowSums(fims_cnaa)
-#   # for (i in 1:length(c(t(em_input$L.age.obs$fleet1)))){
-#   #
-#   #   if (c(t(em_input$L.age.obs$fleet1))[i] == 0) {
-#   #     expect_lte(abs(c(t(fims_cnaa_proportion))[i] - (c(t(em_input$L.age.obs$fleet1))[i]+0.001))/
-#   #                  (c(t(em_input$L.age.obs$fleet1))[i]+0.001), 18)
-#   #   } else {
-#   #     expect_lte(abs(c(t(fims_cnaa_proportion))[i] - c(t(em_input$L.age.obs$fleet1))[i])/
-#   #                  c(t(em_input$L.age.obs$fleet1))[i], 3) # Inf when i = 299; landings was 0.
-#   #   }
-#   # }
-#
-#
-#   # Expected survey index - om_output
-#   index_are <- rep(0, length(om_output$survey_index_biomass$survey1))
-#   for (i in 1:length(om_output$survey_index_biomass$survey1)) {
-#     index_are[i] <- abs(fims_survey[i, 1] - om_output$survey_index_biomass$survey1[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of survey index
-#   expect_lte(
-#     sum(index_are > qnorm(.975) * fims_survey[, 2]),
-#     0.05 * length(om_output$survey_index_biomass$survey1)
-#   )
-#
-#   # Expected survey index - em_input
-#   index_are <- rep(0, length(em_input$surveyB.obs$survey1))
-#   for (i in 1:length(em_input$surveyB.obs$survey1)) {
-#     index_are[i] <- abs(fims_survey[i, 1] - em_input$surveyB.obs$survey1[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of survey index
-#   # expect_lte(
-#   #   sum(index_are > qnorm(.975) * fims_survey[, 2]),
-#   #   0.05 * length(em_input$surveyB.obs$survey1)
-#   # )
-#
-#   for (i in 1:length(em_input$surveyB.obs$survey1)) {
-#     expect_lte(abs(fims_survey[i, 1] - em_input$surveyB.obs$survey1[i]) /
-#       em_input$surveyB.obs$survey1[i], 0.25)
-#   }
-#
-#   # Expected survey number at age
-#   # for (i in 1:length(c(t(om_output$survey_age_comp$survey1)))){
-#   #   expect_lte(abs(report$cnaa[i,2] - c(t(om_output$survey_age_comp$survey1))[i])/
-#   #                c(t(om_output$survey_age_comp$survey1))[i], 0.001)
-#   # }
-#
-#   # Expected catch number at age in proportion
-#   # fims_cnaa <- matrix(report$cnaa[1:(om_input$nyr*om_input$nages), 2],
-#   #                     nrow = om_input$nyr, byrow = TRUE)
-#   # fims_cnaa_proportion <- fims_cnaa/rowSums(fims_cnaa)
-#   #
-#   # for (i in 1:length(c(t(em_input$survey.age.obs)))){
-#   #   expect_lte(abs(c(t(fims_cnaa_proportion))[i] - c(t(em_input$L.age.obs$fleet1))[i])/
-#   #                c(t(em_input$L.age.obs$fleet1))[i], 0.15)
-#   # }
-#
-#
-#   estimation_env$fims$clear()
-# })
+
+test_that("deterministic test of fims", {
+  # run function defined above to set up the test environment
+  deterministic_env <- setup_fims(
+    om_input = om_input,
+    om_output = om_output,
+    em_input = em_input
+  )
+  # Set-up
+  deterministic_env$fims$CreateTMBModel()
+  # CreateTMBModel calls a function in information that loops
+  # over all the populations and fleets and sets all the pointers
+  parameters <- list(p = deterministic_env$fims$get_fixed())
+  # get_fixed function is an Rcpp function that loops over all Rcpp
+  # modules and returned a vector of parameters being estimated
+
+  # Set up TMB's computational graph
+  obj <- MakeADFun(data = list(), parameters, DLL = "FIMS")
+
+  # Calculate standard errors
+  sdr <- TMB::sdreport(obj)
+  sdr_fixed <- summary(sdr, "fixed")
+
+  # Call report using deterministic parameter values
+  # obj$report() requires parameter list to avoid errors
+  report <- obj$report(obj$par)
+
+  # Compare log(R0) to true value
+  fims_logR0 <- sdr_fixed[1, "Estimate"]
+  expect_gt(fims_logR0, 0.0)
+  expect_equal(fims_logR0, log(om_input$R0))
+
+  # Compare numbers at age to true value
+  for (i in 1:length(c(t(om_output$N.age)))) {
+    expect_equal(report$naa[[1]][i], c(t(om_output$N.age))[i])
+  }
+
+  # Compare biomass to true value
+  for (i in 1:length(om_output$biomass.mt)) {
+    expect_equal(report$biomass[[1]][i], om_output$biomass.mt[i])
+  }
+
+  # Compare spawning biomass to true value
+  for (i in 1:length(om_output$SSB)) {
+    expect_equal(report$ssb[[1]][i], om_output$SSB[i])
+  }
+
+  # Compare recruitment to true value
+  fims_naa <- matrix(report$naa[[1]][1:(om_input$nyr * om_input$nages)],
+    nrow = om_input$nyr, byrow = TRUE
+  )
+
+  # loop over years to compare recruitment by year
+  for (i in 1:om_input$nyr) {
+    expect_equal(fims_naa[i, 1], om_output$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$nyr, 1],
+    report$recruitment[[1]][1:om_input$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$nyr) {
+    expect_equal(report$recruitment[[1]][i], om_output$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$logR.resid[-1])
+
+  # F (fixed at initial "true" values)
+  expect_equal(report$F_mort[[1]], om_output$f)
+
+  # Expected catch
+  fims_index <- report$exp_index
+  for (i in 1:length(om_output$L.mt$fleet1)) {
+    expect_equal(fims_index[[1]][i], om_output$L.mt$fleet1[i])
+  }
+
+  # Expect small relative error for deterministic test
+  fims_object_are <- rep(0, length(em_input$L.obs$fleet1))
+  for (i in 1:length(em_input$L.obs$fleet1)) {
+    fims_object_are[i] <- abs(fims_index[[1]][i] - em_input$L.obs$fleet1[i]) / em_input$L.obs$fleet1[i]
+  }
+
+  # Expect 95% of relative error to be within 2*cv
+  expect_lte(sum(fims_object_are > om_input$cv.L$fleet1 * 2.0), length(em_input$L.obs$fleet1) * 0.05)
+
+  # Compare expected catch number at age to true values
+  for (i in 1:length(c(t(om_output$L.age$fleet1)))) {
+    expect_equal(report$cnaa[[1]][i], c(t(om_output$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$nyr * om_input$nages)],
+    nrow = om_input$nyr, byrow = TRUE
+  )
+  fims_cnaa_proportion <- fims_cnaa / rowSums(fims_cnaa)
+  om_cnaa_proportion <- om_output$L.age$fleet1 / rowSums(om_output$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$nyr * om_input$nages)], nrow = om_input$nyr, byrow = TRUE)
+  expect_equal(fims_index[[2]], apply(cwaa, 1, sum) * om_output$survey_q$survey1)
+
+  for (i in 1:length(om_output$survey_index_biomass$survey1)) {
+    expect_equal(fims_index[[2]][i], om_output$survey_index_biomass$survey1[i])
+  }
+
+  fims_object_are <- rep(0, length(em_input$surveyB.obs$survey1))
+  for (i in 1:length(em_input$survey.obs$survey1)) {
+    fims_object_are[i] <- abs(fims_index[[2]][i] - em_input$surveyB.obs$survey1[i]) / em_input$surveyB.obs$survey1[i]
+  }
+  # Expect 95% of relative error to be within 2*cv
+  expect_lte(sum(fims_object_are > om_input$cv.survey$survey1 * 2.0), length(em_input$surveyB.obs$survey1) * 0.05)
+
+  # Expected catch number at age in proportion
+  fims_cnaa <- matrix(report$cnaa[[2]][1:(om_input$nyr * om_input$nages)],
+    nrow = om_input$nyr, byrow = TRUE
+  )
+
+  for (i in 1:length(c(t(om_output$survey_age_comp$survey1)))) {
+    expect_equal(report$cnaa[[2]][i], c(t(om_output$survey_age_comp$survey1))[i])
+  }
+
+  fims_cnaa_proportion <- fims_cnaa / rowSums(fims_cnaa)
+  om_cnaa_proportion <- om_output$survey_age_comp$survey1 / rowSums(om_output$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])
+  }
+  # clear memory
+  deterministic_env$fims$clear()
+})
+
+
+test_that("nll test of fims", {
+  nll_env <- setup_fims(
+    om_input = om_input,
+    om_output = om_output,
+    em_input = em_input
+  )
+  # Set-up TMB
+  nll_env$fims$CreateTMBModel()
+  parameters <- list(p = nll_env$fims$get_fixed())
+  par_list <- 1:length(parameters[[1]])
+  par_list[2:length(par_list)] <- NA
+  map <- list(p = factor(par_list))
+
+  obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS", map = map)
+
+  sdr <- TMB::sdreport(obj)
+  sdr_fixed <- summary(sdr, "fixed")
+
+  # log(R0)
+  fims_logR0 <- sdr_fixed[1, "Estimate"]
+  # expect_lte(abs(fims_logR0 - log(om_input$R0)) / log(om_input$R0), 0.0001)
+  expect_equal(fims_logR0, log(om_input$R0))
+
+  # Call report using deterministic parameter values
+  # obj$report() requires parameter list to avoid errors
+  report <- obj$report(obj$par)
+  obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS", map = map)
+  jnll <- obj$fn()
+
+
+  # recruitment likelihood
+  # log_devs is of length nyr-1
+  rec_nll <- 0
+  for(i in 1:(om_input$nyr - 1)){
+    rec_nll <- rec_nll - dnorm(
+      nll_env$recruitment$log_devs[i]$value, 0,
+      om_input$logR_sd, TRUE
+    )
+  }
+
+  # catch and survey index expected likelihoods
+  index_nll_fleet <- -sum(dnorm(
+    log(nll_env$catch),
+    log(om_output$L.mt$fleet1),
+    sqrt(log(em_input$cv.L$fleet1^2 + 1)), TRUE
+  ))
+  index_nll_survey <- -sum(dnorm(
+    log(nll_env$survey_index),
+    log(om_output$survey_index_biomass$survey1),
+    sqrt(log(em_input$cv.survey$survey1^2 + 1)), TRUE
+  ))
+  index_nll <- index_nll_fleet + index_nll_survey
+  # age comp likelihoods
+  fishing_acomp_observed <- em_input$L.age.obs$fleet1
+  fishing_acomp_expected <- om_output$L.age$fleet1 / rowSums(om_output$L.age$fleet1)
+  survey_acomp_observed <- em_input$survey.age.obs$survey1
+  survey_acomp_expected <- om_output$survey_age_comp$survey1 / rowSums(om_output$survey_age_comp$survey1)
+  age_comp_nll_fleet <- age_comp_nll_survey <- 0
+  for (y in 1:om_input$nyr) {
+    age_comp_nll_fleet <- age_comp_nll_fleet -
+      dmultinom(
+        fishing_acomp_observed[y, ] * em_input$n.L$fleet1, em_input$n.L$fleet1,
+        fishing_acomp_expected[y, ], TRUE
+      )
+
+    age_comp_nll_survey <- age_comp_nll_survey -
+      dmultinom(
+        survey_acomp_observed[y, ] * em_input$n.survey$survey1, em_input$n.survey$survey1,
+        survey_acomp_expected[y, ], TRUE
+      )
+  }
+  age_comp_nll <- age_comp_nll_fleet + age_comp_nll_survey
+  expected_jnll <- rec_nll + index_nll + age_comp_nll
+
+  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)
+
+  nll_env$fims$clear()
+})
+
+test_that("estimation test of fims", {
+  estimation_env <- setup_fims(
+    om_input = om_input,
+    om_output = om_output,
+    em_input = em_input
+  )
+
+  # Set-up TMB
+  estimation_env$fims$CreateTMBModel()
+  # Create parameter list from Rcpp modules
+  parameters <- list(p = estimation_env$fims$get_fixed())
+  obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS")
+
+  opt <- with(obj, optim(par, fn, gr,
+    method = "BFGS",
+    control = list(maxit = 1000000, reltol = 1e-15)
+  ))
+
+  # Call report using MLE parameter values
+  # obj$report() requires parameter list to avoid errors
+  report <- obj$report(obj$env$last.par.best)
+  sdr <- TMB::sdreport(obj)
+  sdr_report <- summary(sdr, "report")
+  # Numbers at age
+  # Estimates and SE for NAA
+  sdr_naa <- sdr_report[which(rownames(sdr_report) == "NAA"), ]
+  naa_are <- rep(0, length(c(t(om_output$N.age))))
+  for (i in 1:length(c(t(om_output$N.age)))) {
+    naa_are[i] <- abs(sdr_naa[i, 1] - c(t(om_output$N.age))[i])
+  }
+  # Expect 95% of absolute error to be within 2*SE of NAA
+  expect_lte(
+    sum(naa_are > qnorm(.975) * sdr_naa[1:length(c(t(om_output$N.age))), 2]),
+    0.05 * length(c(t(om_output$N.age)))
+  )
+
+  # Biomass
+  sdr_biomass <- sdr_report[which(rownames(sdr_report) == "Biomass"), ]
+  biomass_are <- rep(0, length(om_output$biomass.mt))
+  for (i in 1:length(om_output$biomass.mt)) {
+    biomass_are[i] <- abs(sdr_biomass[i, 1] - om_output$biomass.mt[i]) # / om_output$biomass.mt[i]
+    # expect_lte(biomass_are[i], 0.15)
+  }
+  expect_lte(
+    sum(biomass_are > qnorm(.975) * sdr_biomass[1:length(om_output$biomass.mt), 2]),
+    0.05 * length(om_output$biomass.mt)
+  )
+
+  # Spawning biomass
+  sdr_sb <- sdr_report[which(rownames(sdr_report) == "SSB"), ]
+  sb_are <- rep(0, length(om_output$SSB))
+  for (i in 1:length(om_output$SSB)) {
+    sb_are[i] <- abs(sdr_sb[i, 1] - om_output$SSB[i]) # / om_output$SSB[i]
+    # expect_lte(sb_are[i], 0.15)
+  }
+  expect_lte(
+    sum(sb_are > qnorm(.975) * sdr_sb[1:length(om_output$SSB), 2]),
+    0.05 * length(om_output$SSB)
+  )
+
+  # Recruitment
+  fims_naa <- matrix(report$naa[[1]][1:(om_input$nyr * om_input$nages)],
+    nrow = om_input$nyr, byrow = TRUE
+  )
+  sdr_naa1_vec <- sdr_report[which(rownames(sdr_report) == "NAA"), 2]
+  sdr_naa1 <- sdr_naa1_vec[seq(1, om_input$nyr * om_input$nages, by = om_input$nages)]
+  fims_naa1_are <- rep(0, om_input$nyr)
+  for (i in 1:om_input$nyr) {
+    fims_naa1_are[i] <- abs(fims_naa[i, 1] - om_output$N.age[i, 1]) # /
+    # om_output$N.age[i, 1]
+    # expect_lte(fims_naa1_are[i], 0.25)
+  }
+  expect_lte(
+    sum(fims_naa1_are > qnorm(.975) * sdr_naa1[1:length(om_output$SSB)]),
+    0.05 * length(om_output$SSB)
+  )
+
+  expect_equal(
+    fims_naa[, 1],
+    report$recruitment[[1]][1:om_input$nyr]
+  )
+
+  # recruitment log deviations
+  # the initial value of om_input$logR.resid is dropped from the model
+  sdr_rdev <- sdr_report[which(rownames(sdr_report) == "LogRecDev"), ]
+  rdev_are <- rep(0, length(om_input$logR.resid) - 1)
+
+  for (i in 1:(length(report$log_recruit_dev[[1]]) - 1)) {
+    rdev_are[i] <- abs(report$log_recruit_dev[[1]][i] - om_input$logR.resid[i + 1]) # /
+    #   exp(om_input$logR.resid[i])
+    # expect_lte(rdev_are[i], 1) # 1
+  }
+  expect_lte(
+    sum(rdev_are > qnorm(.975) * sdr_rdev[1:length(om_input$logR.resid) - 1, 2]),
+    0.05 * length(om_input$logR.resid)
+  )
+
+  # F (needs to be updated when std.error is available)
+  sdr_F <- sdr_report[which(rownames(sdr_report) == "FMort"), ]
+  f_are <- rep(0, length(om_output$f))
+  for (i in 1:length(om_output$f)) {
+    f_are[i] <- abs(sdr_F[i, 1] - om_output$f[i])
+  }
+  # Expect 95% of absolute error to be within 2*SE of Fmort
+  expect_lte(
+    sum(f_are > qnorm(.975) * sdr_F[1:length(om_output$f), 2]),
+    0.05 * length(om_output$f)
+  )
+
+  # Expected fishery catch and survey index
+  fims_index <- sdr_report[which(rownames(sdr_report) == "ExpectedIndex"), ]
+  fims_catch <- fims_index[1:om_input$nyr, ]
+  fims_survey <- fims_index[(om_input$nyr + 1):(om_input$nyr * 2), ]
+
+  # Expected fishery catch - om_output
+  catch_are <- rep(0, length(om_output$L.mt$fleet1))
+  for (i in 1:length(om_output$L.mt$fleet1)) {
+    catch_are[i] <- abs(fims_catch[i, 1] - om_output$L.mt$fleet1[i])
+  }
+  # Expect 95% of absolute error to be within 2*SE of fishery catch
+  expect_lte(
+    sum(catch_are > qnorm(.975) * fims_catch[, 2]),
+    0.05 * length(om_output$L.mt$fleet1)
+  )
+
+  # Expected fishery catch - em_input
+  catch_are <- rep(0, length(em_input$L.obs$fleet1))
+  for (i in 1:length(em_input$L.obs$fleet1)) {
+    catch_are[i] <- abs(fims_catch[i, 1] - em_input$L.obs$fleet1[i])
+  }
+  # Expect 95% of absolute error to be within 2*SE of fishery catch
+  expect_lte(
+    sum(catch_are > qnorm(.975) * fims_catch[, 2]),
+    0.05 * length(em_input$L.obs$fleet1)
+  )
+
+
+  # Expected fishery catch number at age
+  sdr_cnaa <- sdr_report[which(rownames(sdr_report) == "CNAA"), ]
+  cnaa_are <- rep(0, length(c(t(om_output$L.age$fleet1))))
+  for (i in 1:length(c(t(om_output$L.age$fleet1)))) {
+    cnaa_are[i] <- abs(sdr_cnaa[i, 1] - c(t(om_output$L.age$fleet1))[i])
+  }
+  # Expect 95% of absolute error to be within 2*SE of CNAA
+  expect_lte(
+    sum(cnaa_are > qnorm(.975) * sdr_cnaa[, 2]),
+    0.05 * length(c(t(om_output$L.age$fleet1)))
+  )
+
+
+  # Expected catch number at age in proportion
+  # fims_cnaa <- matrix(report$cnaa[1:(om_input$nyr*om_input$nages), 1],
+  #                     nrow = om_input$nyr, byrow = TRUE)
+  # fims_cnaa_proportion <- fims_cnaa/rowSums(fims_cnaa)
+  # for (i in 1:length(c(t(em_input$L.age.obs$fleet1)))){
+  #
+  #   if (c(t(em_input$L.age.obs$fleet1))[i] == 0) {
+  #     expect_lte(abs(c(t(fims_cnaa_proportion))[i] - (c(t(em_input$L.age.obs$fleet1))[i]+0.001))/
+  #                  (c(t(em_input$L.age.obs$fleet1))[i]+0.001), 18)
+  #   } else {
+  #     expect_lte(abs(c(t(fims_cnaa_proportion))[i] - c(t(em_input$L.age.obs$fleet1))[i])/
+  #                  c(t(em_input$L.age.obs$fleet1))[i], 3) # Inf when i = 299; landings was 0.
+  #   }
+  # }
+
+
+  # Expected survey index - om_output
+  index_are <- rep(0, length(om_output$survey_index_biomass$survey1))
+  for (i in 1:length(om_output$survey_index_biomass$survey1)) {
+    index_are[i] <- abs(fims_survey[i, 1] - om_output$survey_index_biomass$survey1[i])
+  }
+  # Expect 95% of absolute error to be within 2*SE of survey index
+  expect_lte(
+    sum(index_are > qnorm(.975) * fims_survey[, 2]),
+    0.05 * length(om_output$survey_index_biomass$survey1)
+  )
+
+  # Expected survey index - em_input
+  index_are <- rep(0, length(em_input$surveyB.obs$survey1))
+  for (i in 1:length(em_input$surveyB.obs$survey1)) {
+    index_are[i] <- abs(fims_survey[i, 1] - em_input$surveyB.obs$survey1[i])
+  }
+  # Expect 95% of absolute error to be within 2*SE of survey index
+  # expect_lte(
+  #   sum(index_are > qnorm(.975) * fims_survey[, 2]),
+  #   0.05 * length(em_input$surveyB.obs$survey1)
+  # )
+
+  for (i in 1:length(em_input$surveyB.obs$survey1)) {
+    expect_lte(abs(fims_survey[i, 1] - em_input$surveyB.obs$survey1[i]) /
+      em_input$surveyB.obs$survey1[i], 0.25)
+  }
+
+  # Expected survey number at age
+  # for (i in 1:length(c(t(om_output$survey_age_comp$survey1)))){
+  #   expect_lte(abs(report$cnaa[i,2] - c(t(om_output$survey_age_comp$survey1))[i])/
+  #                c(t(om_output$survey_age_comp$survey1))[i], 0.001)
+  # }
+
+  # Expected catch number at age in proportion
+  # fims_cnaa <- matrix(report$cnaa[1:(om_input$nyr*om_input$nages), 2],
+  #                     nrow = om_input$nyr, byrow = TRUE)
+  # fims_cnaa_proportion <- fims_cnaa/rowSums(fims_cnaa)
+  #
+  # for (i in 1:length(c(t(em_input$survey.age.obs)))){
+  #   expect_lte(abs(c(t(fims_cnaa_proportion))[i] - c(t(em_input$L.age.obs$fleet1))[i])/
+  #                c(t(em_input$L.age.obs$fleet1))[i], 0.15)
+  # }
+
+
+  estimation_env$fims$clear()
+})
 
 test_that("run FIMS in a for loop with missing values", {
   for (i in 1:5) {
@@ -820,189 +825,189 @@ test_that("run FIMS in a for loop with missing values", {
   }
 })
 
-# test_that("agecomp in proportion works", {
-#   n.L_original <- om_input$n.L$fleet1
-#   n.survey_original <- om_input$n.survey$survey1
-#   om_input$n.L$fleet1 <- 1
-#   om_input$n.survey$survey1 <- 1
-#   on.exit(om_input$n.L$fleet1 <- n.L_original, add = TRUE)
-#   on.exit(om_input$n.survey$survey1 <- n.survey_original, add = TRUE)
-#
-#   agecomp_in_proportion_env <- setup_fims(
-#     om_input = om_input,
-#     om_output = om_output,
-#     em_input = em_input
-#   )
-#
-#   # Set-up TMB
-#   agecomp_in_proportion_env$fims$CreateTMBModel()
-#   # Create parameter list from Rcpp modules
-#   parameters <- list(p = agecomp_in_proportion_env$fims$get_fixed())
-#   obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS")
-#
-#   opt <- with(obj, optim(par, fn, gr,
-#     method = "BFGS",
-#     control = list(maxit = 1000000, reltol = 1e-15)
-#   ))
-#
-#   # Call report using MLE parameter values
-#   # obj$report() requires parameter list to avoid errors
-#   report <- obj$report(obj$env$last.par.best)
-#   sdr <- TMB::sdreport(obj)
-#   sdr_report <- summary(sdr, "report")
-#   # Numbers at age
-#   # Estimates and SE for NAA
-#   sdr_naa <- sdr_report[which(rownames(sdr_report) == "NAA"), ]
-#   naa_are <- rep(0, length(c(t(om_output$N.age))))
-#   for (i in 1:length(c(t(om_output$N.age)))) {
-#     naa_are[i] <- abs(sdr_naa[i, 1] - c(t(om_output$N.age))[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of NAA
-#   expect_lte(
-#     sum(naa_are > qnorm(.975) * sdr_naa[1:length(c(t(om_output$N.age))), 2]),
-#     0.05 * length(c(t(om_output$N.age)))
-#   )
-#
-#   # Biomass
-#   sdr_biomass <- sdr_report[which(rownames(sdr_report) == "Biomass"), ]
-#   biomass_are <- rep(0, length(om_output$biomass.mt))
-#   for (i in 1:length(om_output$biomass.mt)) {
-#     biomass_are[i] <- abs(sdr_biomass[i, 1] - om_output$biomass.mt[i]) # / om_output$biomass.mt[i]
-#     # expect_lte(biomass_are[i], 0.15)
-#   }
-#   expect_lte(
-#     sum(biomass_are > qnorm(.975) * sdr_biomass[1:length(om_output$biomass.mt), 2]),
-#     0.05 * length(om_output$biomass.mt)
-#   )
-#
-#   # Spawning biomass
-#   sdr_sb <- sdr_report[which(rownames(sdr_report) == "SSB"), ]
-#   sb_are <- rep(0, length(om_output$SSB))
-#   for (i in 1:length(om_output$SSB)) {
-#     sb_are[i] <- abs(sdr_sb[i, 1] - om_output$SSB[i]) # / om_output$SSB[i]
-#     # expect_lte(sb_are[i], 0.15)
-#   }
-#   expect_lte(
-#     sum(sb_are > qnorm(.975) * sdr_sb[1:length(om_output$SSB), 2]),
-#     0.05 * length(om_output$SSB)
-#   )
-#
-#   # Recruitment
-#   fims_naa <- matrix(report$naa[[1]][1:(om_input$nyr * om_input$nages)],
-#     nrow = om_input$nyr, byrow = TRUE
-#   )
-#   sdr_naa1_vec <- sdr_report[which(rownames(sdr_report) == "NAA"), 2]
-#   sdr_naa1 <- sdr_naa1_vec[seq(1, om_input$nyr * om_input$nages, by = om_input$nages)]
-#   fims_naa1_are <- rep(0, om_input$nyr)
-#   for (i in 1:om_input$nyr) {
-#     fims_naa1_are[i] <- abs(fims_naa[i, 1] - om_output$N.age[i, 1]) # /
-#     # om_output$N.age[i, 1]
-#     # expect_lte(fims_naa1_are[i], 0.25)
-#   }
-#   expect_lte(
-#     sum(fims_naa1_are > qnorm(.975) * sdr_naa1[1:length(om_output$SSB)]),
-#     0.05 * length(om_output$SSB)
-#   )
-#
-#   expect_equal(
-#     fims_naa[, 1],
-#     report$recruitment[[1]][1:om_input$nyr]
-#   )
-#
-#   # recruitment log deviations
-#   # the initial value of om_input$logR.resid is dropped from the model
-#   sdr_rdev <- sdr_report[which(rownames(sdr_report) == "LogRecDev"), ]
-#   rdev_are <- rep(0, length(om_input$logR.resid) - 1)
-#
-#   for (i in 1:(length(report$log_recruit_dev[[1]]) - 1)) {
-#     rdev_are[i] <- abs(report$log_recruit_dev[[1]][i] - om_input$logR.resid[i + 1]) # /
-#     #   exp(om_input$logR.resid[i])
-#     # expect_lte(rdev_are[i], 1) # 1
-#   }
-#   expect_lte(
-#     sum(rdev_are > qnorm(.975) * sdr_rdev[1:length(om_input$logR.resid) - 1, 2]),
-#     0.05 * length(om_input$logR.resid)
-#   )
-#
-#   # F (needs to be updated when std.error is available)
-#   sdr_F <- sdr_report[which(rownames(sdr_report) == "FMort"), ]
-#   f_are <- rep(0, length(om_output$f))
-#   for (i in 1:length(om_output$f)) {
-#     f_are[i] <- abs(sdr_F[i, 1] - om_output$f[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of Fmort
-#   expect_lte(
-#     sum(f_are > qnorm(.975) * sdr_F[1:length(om_output$f), 2]),
-#     0.05 * length(om_output$f)
-#   )
-#
-#   # Expected fishery catch and survey index
-#   fims_index <- sdr_report[which(rownames(sdr_report) == "ExpectedIndex"), ]
-#   fims_catch <- fims_index[1:om_input$nyr, ]
-#   fims_survey <- fims_index[(om_input$nyr + 1):(om_input$nyr * 2), ]
-#
-#   # Expected fishery catch - om_output
-#   catch_are <- rep(0, length(om_output$L.mt$fleet1))
-#   for (i in 1:length(om_output$L.mt$fleet1)) {
-#     catch_are[i] <- abs(fims_catch[i, 1] - om_output$L.mt$fleet1[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of fishery catch
-#   expect_lte(
-#     sum(catch_are > qnorm(.975) * fims_catch[, 2]),
-#     0.05 * length(om_output$L.mt$fleet1)
-#   )
-#
-#   # Expected fishery catch - em_input
-#   catch_are <- rep(0, length(em_input$L.obs$fleet1))
-#   for (i in 1:length(em_input$L.obs$fleet1)) {
-#     catch_are[i] <- abs(fims_catch[i, 1] - em_input$L.obs$fleet1[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of fishery catch
-#   expect_lte(
-#     sum(catch_are > qnorm(.975) * fims_catch[, 2]),
-#     0.05 * length(em_input$L.obs$fleet1)
-#   )
-#
-#
-#   # Expected fishery catch number at age
-#   sdr_cnaa <- sdr_report[which(rownames(sdr_report) == "CNAA"), ]
-#   cnaa_are <- rep(0, length(c(t(om_output$L.age$fleet1))))
-#   for (i in 1:length(c(t(om_output$L.age$fleet1)))) {
-#     cnaa_are[i] <- abs(sdr_cnaa[i, 1] - c(t(om_output$L.age$fleet1))[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of CNAA
-#   expect_lte(
-#     sum(cnaa_are > qnorm(.975) * sdr_cnaa[, 2]),
-#     0.05 * length(c(t(om_output$L.age$fleet1)))
-#   )
-#   print(cbind(sdr_cnaa, c(t(om_output$L.age$fleet1))))
-#
-#   # Expected survey index - om_output
-#   index_are <- rep(0, length(om_output$survey_index_biomass$survey1))
-#   for (i in 1:length(om_output$survey_index_biomass$survey1)) {
-#     index_are[i] <- abs(fims_survey[i, 1] - om_output$survey_index_biomass$survey1[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of survey index
-#   expect_lte(
-#     sum(index_are > qnorm(.975) * fims_survey[, 2]),
-#     0.05 * length(om_output$survey_index_biomass$survey1)
-#   )
-#
-#   # Expected survey index - em_input
-#   index_are <- rep(0, length(em_input$surveyB.obs$survey1))
-#   for (i in 1:length(em_input$surveyB.obs$survey1)) {
-#     index_are[i] <- abs(fims_survey[i, 1] - em_input$surveyB.obs$survey1[i])
-#   }
-#   # Expect 95% of absolute error to be within 2*SE of survey index
-#   # expect_lte(
-#   #   sum(index_are > qnorm(.975) * fims_survey[, 2]),
-#   #   0.05 * length(em_input$surveyB.obs$survey1)
-#   # )
-#
-#   for (i in 1:length(em_input$surveyB.obs$survey1)) {
-#     expect_lte(abs(fims_survey[i, 1] - em_input$surveyB.obs$survey1[i]) /
-#       em_input$surveyB.obs$survey1[i], 0.25)
-#   }
-#
-#   agecomp_in_proportion_env$fims$clear()
-# })
+test_that("agecomp in proportion works", {
+  n.L_original <- om_input$n.L$fleet1
+  n.survey_original <- om_input$n.survey$survey1
+  om_input$n.L$fleet1 <- 1
+  om_input$n.survey$survey1 <- 1
+  on.exit(om_input$n.L$fleet1 <- n.L_original, add = TRUE)
+  on.exit(om_input$n.survey$survey1 <- n.survey_original, add = TRUE)
+
+  agecomp_in_proportion_env <- setup_fims(
+    om_input = om_input,
+    om_output = om_output,
+    em_input = em_input
+  )
+
+  # Set-up TMB
+  agecomp_in_proportion_env$fims$CreateTMBModel()
+  # Create parameter list from Rcpp modules
+  parameters <- list(p = agecomp_in_proportion_env$fims$get_fixed())
+  obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS")
+
+  opt <- with(obj, optim(par, fn, gr,
+    method = "BFGS",
+    control = list(maxit = 1000000, reltol = 1e-15)
+  ))
+
+  # Call report using MLE parameter values
+  # obj$report() requires parameter list to avoid errors
+  report <- obj$report(obj$env$last.par.best)
+  sdr <- TMB::sdreport(obj)
+  sdr_report <- summary(sdr, "report")
+  # Numbers at age
+  # Estimates and SE for NAA
+  sdr_naa <- sdr_report[which(rownames(sdr_report) == "NAA"), ]
+  naa_are <- rep(0, length(c(t(om_output$N.age))))
+  for (i in 1:length(c(t(om_output$N.age)))) {
+    naa_are[i] <- abs(sdr_naa[i, 1] - c(t(om_output$N.age))[i])
+  }
+  # Expect 95% of absolute error to be within 2*SE of NAA
+  expect_lte(
+    sum(naa_are > qnorm(.975) * sdr_naa[1:length(c(t(om_output$N.age))), 2]),
+    0.05 * length(c(t(om_output$N.age)))
+  )
+
+  # Biomass
+  sdr_biomass <- sdr_report[which(rownames(sdr_report) == "Biomass"), ]
+  biomass_are <- rep(0, length(om_output$biomass.mt))
+  for (i in 1:length(om_output$biomass.mt)) {
+    biomass_are[i] <- abs(sdr_biomass[i, 1] - om_output$biomass.mt[i]) # / om_output$biomass.mt[i]
+    # expect_lte(biomass_are[i], 0.15)
+  }
+  expect_lte(
+    sum(biomass_are > qnorm(.975) * sdr_biomass[1:length(om_output$biomass.mt), 2]),
+    0.05 * length(om_output$biomass.mt)
+  )
+
+  # Spawning biomass
+  sdr_sb <- sdr_report[which(rownames(sdr_report) == "SSB"), ]
+  sb_are <- rep(0, length(om_output$SSB))
+  for (i in 1:length(om_output$SSB)) {
+    sb_are[i] <- abs(sdr_sb[i, 1] - om_output$SSB[i]) # / om_output$SSB[i]
+    # expect_lte(sb_are[i], 0.15)
+  }
+  expect_lte(
+    sum(sb_are > qnorm(.975) * sdr_sb[1:length(om_output$SSB), 2]),
+    0.05 * length(om_output$SSB)
+  )
+
+  # Recruitment
+  fims_naa <- matrix(report$naa[[1]][1:(om_input$nyr * om_input$nages)],
+    nrow = om_input$nyr, byrow = TRUE
+  )
+  sdr_naa1_vec <- sdr_report[which(rownames(sdr_report) == "NAA"), 2]
+  sdr_naa1 <- sdr_naa1_vec[seq(1, om_input$nyr * om_input$nages, by = om_input$nages)]
+  fims_naa1_are <- rep(0, om_input$nyr)
+  for (i in 1:om_input$nyr) {
+    fims_naa1_are[i] <- abs(fims_naa[i, 1] - om_output$N.age[i, 1]) # /
+    # om_output$N.age[i, 1]
+    # expect_lte(fims_naa1_are[i], 0.25)
+  }
+  expect_lte(
+    sum(fims_naa1_are > qnorm(.975) * sdr_naa1[1:length(om_output$SSB)]),
+    0.05 * length(om_output$SSB)
+  )
+
+  expect_equal(
+    fims_naa[, 1],
+    report$recruitment[[1]][1:om_input$nyr]
+  )
+
+  # recruitment log deviations
+  # the initial value of om_input$logR.resid is dropped from the model
+  sdr_rdev <- sdr_report[which(rownames(sdr_report) == "LogRecDev"), ]
+  rdev_are <- rep(0, length(om_input$logR.resid) - 1)
+
+  for (i in 1:(length(report$log_recruit_dev[[1]]) - 1)) {
+    rdev_are[i] <- abs(report$log_recruit_dev[[1]][i] - om_input$logR.resid[i + 1]) # /
+    #   exp(om_input$logR.resid[i])
+    # expect_lte(rdev_are[i], 1) # 1
+  }
+  expect_lte(
+    sum(rdev_are > qnorm(.975) * sdr_rdev[1:length(om_input$logR.resid) - 1, 2]),
+    0.05 * length(om_input$logR.resid)
+  )
+
+  # F (needs to be updated when std.error is available)
+  sdr_F <- sdr_report[which(rownames(sdr_report) == "FMort"), ]
+  f_are <- rep(0, length(om_output$f))
+  for (i in 1:length(om_output$f)) {
+    f_are[i] <- abs(sdr_F[i, 1] - om_output$f[i])
+  }
+  # Expect 95% of absolute error to be within 2*SE of Fmort
+  expect_lte(
+    sum(f_are > qnorm(.975) * sdr_F[1:length(om_output$f), 2]),
+    0.05 * length(om_output$f)
+  )
+
+  # Expected fishery catch and survey index
+  fims_index <- sdr_report[which(rownames(sdr_report) == "ExpectedIndex"), ]
+  fims_catch <- fims_index[1:om_input$nyr, ]
+  fims_survey <- fims_index[(om_input$nyr + 1):(om_input$nyr * 2), ]
+
+  # Expected fishery catch - om_output
+  catch_are <- rep(0, length(om_output$L.mt$fleet1))
+  for (i in 1:length(om_output$L.mt$fleet1)) {
+    catch_are[i] <- abs(fims_catch[i, 1] - om_output$L.mt$fleet1[i])
+  }
+  # Expect 95% of absolute error to be within 2*SE of fishery catch
+  expect_lte(
+    sum(catch_are > qnorm(.975) * fims_catch[, 2]),
+    0.05 * length(om_output$L.mt$fleet1)
+  )
+
+  # Expected fishery catch - em_input
+  catch_are <- rep(0, length(em_input$L.obs$fleet1))
+  for (i in 1:length(em_input$L.obs$fleet1)) {
+    catch_are[i] <- abs(fims_catch[i, 1] - em_input$L.obs$fleet1[i])
+  }
+  # Expect 95% of absolute error to be within 2*SE of fishery catch
+  expect_lte(
+    sum(catch_are > qnorm(.975) * fims_catch[, 2]),
+    0.05 * length(em_input$L.obs$fleet1)
+  )
+
+
+  # Expected fishery catch number at age
+  sdr_cnaa <- sdr_report[which(rownames(sdr_report) == "CNAA"), ]
+  cnaa_are <- rep(0, length(c(t(om_output$L.age$fleet1))))
+  for (i in 1:length(c(t(om_output$L.age$fleet1)))) {
+    cnaa_are[i] <- abs(sdr_cnaa[i, 1] - c(t(om_output$L.age$fleet1))[i])
+  }
+  # Expect 95% of absolute error to be within 2*SE of CNAA
+  expect_lte(
+    sum(cnaa_are > qnorm(.975) * sdr_cnaa[, 2]),
+    0.05 * length(c(t(om_output$L.age$fleet1)))
+  )
+  print(cbind(sdr_cnaa, c(t(om_output$L.age$fleet1))))
+
+  # Expected survey index - om_output
+  index_are <- rep(0, length(om_output$survey_index_biomass$survey1))
+  for (i in 1:length(om_output$survey_index_biomass$survey1)) {
+    index_are[i] <- abs(fims_survey[i, 1] - om_output$survey_index_biomass$survey1[i])
+  }
+  # Expect 95% of absolute error to be within 2*SE of survey index
+  expect_lte(
+    sum(index_are > qnorm(.975) * fims_survey[, 2]),
+    0.05 * length(om_output$survey_index_biomass$survey1)
+  )
+
+  # Expected survey index - em_input
+  index_are <- rep(0, length(em_input$surveyB.obs$survey1))
+  for (i in 1:length(em_input$surveyB.obs$survey1)) {
+    index_are[i] <- abs(fims_survey[i, 1] - em_input$surveyB.obs$survey1[i])
+  }
+  # Expect 95% of absolute error to be within 2*SE of survey index
+  # expect_lte(
+  #   sum(index_are > qnorm(.975) * fims_survey[, 2]),
+  #   0.05 * length(em_input$surveyB.obs$survey1)
+  # )
+
+  for (i in 1:length(em_input$surveyB.obs$survey1)) {
+    expect_lte(abs(fims_survey[i, 1] - em_input$surveyB.obs$survey1[i]) /
+      em_input$surveyB.obs$survey1[i], 0.25)
+  }
+
+  agecomp_in_proportion_env$fims$clear()
+})
diff --git a/tests/testthat/test-rcpp-tmb-distributions.R b/tests/testthat/test-rcpp-tmb-distributions.R
index 1fcf83729..1dd62896d 100644
--- a/tests/testthat/test-rcpp-tmb-distributions.R
+++ b/tests/testthat/test-rcpp-tmb-distributions.R
@@ -15,7 +15,6 @@ test_that("normal_lpdf", {
   dnorm_$x <- new(ParameterVector, y, 1)
   dnorm_$expected_values <- new(ParameterVector, 0, 1)
   dnorm_$log_sd <- new(ParameterVector, log(1), 1)
-  dnorm_$is_na <- FALSE
   # evaluate the density and compare with R
   expect_equal(dnorm_$evaluate(), stats::dnorm(y, 0, 1, TRUE))
   clear()
@@ -31,7 +30,6 @@ test_that("normal_lpdf", {
   dnorm_$x <- new(ParameterVector, y, 10)
   dnorm_$expected_values <- new(ParameterVector, 0, 1)
   dnorm_$log_sd <- new(ParameterVector, log(1), 1)
-  dnorm_$is_na <- rep(FALSE, 10)
   # evaluate the density and compare with R
   expect_equal(dnorm_$evaluate(), sum(stats::dnorm(y, 0, 1, TRUE)))
   clear()
@@ -45,9 +43,8 @@ test_that("normal_lpdf", {
   dnorm_ <- new(TMBDnormDistribution)
   # populate class members
   dnorm_$x <- new(ParameterVector, y, 10)
-  dnorm_$expected_values <- new(ParameterVector, 0.0, 10)
-  dnorm_$log_sd <- new(ParameterVector, log(1), 10)
-  dnorm_$is_na <- rep(FALSE, 10)
+  dnorm_$expected_values <- new(ParameterVector, rep(0,10), 10)
+  dnorm_$log_sd <- new(ParameterVector, rep(log(1),10), 10)
   # evaluate the density and compare with R
   expect_equal(dnorm_$evaluate(), sum(stats::dnorm(y, 0, 1, TRUE)))
   clear()
@@ -62,7 +59,6 @@ test_that("normal_lpdf", {
   # dnorm_$x <- new(FIMS:::ParameterVector, y, 10)
   # dnorm_$expected_values <- new(FIMS:::ParameterVector, 0, 11)
   # dnorm_$log_sd <- new(FIMS:::ParameterVector, log(1), 3)
-  # dnorm_$is_na <- rep(FALSE, 10)
   # clear()
 })
 
@@ -83,10 +79,8 @@ test_that("lognormal_lpdf", {
   dlnorm_$x <- new(ParameterVector, y, 1)
   dlnorm_$expected_values <- new(ParameterVector, 0, 1)
   dlnorm_$log_logsd <- new(ParameterVector, log(1), 1)
-  dlnorm_$is_na <- FALSE
-  dlnorm_$input_type <- "data"
   # evaluate the density and compare with R
-  expect_equal(dlnorm_$evaluate(), stats::dlnorm(y, 0, 1, TRUE))
+  expect_equal(dlnorm_$evaluate(), stats::dlnorm(y, 0, 1, TRUE) + log(y))
   clear()
 
   ## A vector of state variables, but scalar arguments, e.g., a
@@ -100,10 +94,8 @@ test_that("lognormal_lpdf", {
   dlnorm_$x <- new(ParameterVector, y, 10)
   dlnorm_$expected_values <- new(ParameterVector, 0, 1)
   dlnorm_$log_logsd <- new(ParameterVector, log(1), 1)
-  dlnorm_$is_na <- rep(FALSE, 10)
-  dlnorm_$input_type <- "data"
   # evaluate the density and compare with R
-  expect_equal(dlnorm_$evaluate(), sum(stats::dlnorm(y, 0, 1, TRUE)))
+  expect_equal(dlnorm_$evaluate(), sum(stats::dlnorm(y, 0, 1, TRUE)) + sum(log(y)))
   clear()
 
 
@@ -116,12 +108,10 @@ test_that("lognormal_lpdf", {
   dlnorm_ <- new(TMBDlnormDistribution)
   # populate class members
   dlnorm_$x <- new(ParameterVector, y, 10)
-  dlnorm_$expected_values <- new(ParameterVector, 0, 10)
-  dlnorm_$log_logsd <- new(ParameterVector, log(1), 10)
-  dlnorm_$is_na <- rep(FALSE, 10)
-  dlnorm_$input_type <- "data"
+  dlnorm_$expected_values <- new(ParameterVector, rep(0,10), 10)
+  dlnorm_$log_logsd <- new(ParameterVector, rep(log(1),10), 10)
   # evaluate the density and compare with R
-  expect_equal(dlnorm_$evaluate(), sum(stats::dlnorm(y, 0, 1, TRUE)))
+  expect_equal(dlnorm_$evaluate(), sum(stats::dlnorm(y, 0, 1, TRUE)) + sum(log(y)))
   clear()
 
   ## It should error out when there is a dimension mismatch
@@ -135,67 +125,32 @@ test_that("lognormal_lpdf", {
   # dlnorm_$x <- new(ParameterVector, y, 10)
   # dlnorm_$expected_values <- new(ParameterVector, 0, 11)
   # dlnorm_$log_logsd <- new(ParameterVector, log(1), 3)
-  # dlnorm_$is_na <- rep(FALSE, 10)
   # clear()
 
 
 })
+# test not working
+# test_that("multinomial_lpdf", {
+#   # generate data using R stats:rnorm
+#   set.seed(123)
+#   p <- (1:10) / sum(1:10)
+#   x <- stats::rmultinom(1, 100, p)
+#   compdata <- methods::new(AgeComp, 1, 10)
+#   compdata$age_comp_data <- x
+#   # create a fims Rcpp object
+#   # initialize the Dmultinom module
+#   dmultinom_ <- new(TMBDmultinomDistribution)
+#   # populate class members
+#   dmultinom_$expected_values <- new(ParameterVector, p, 10)
+#   dmultinom_$set_observed_data(compdata$get_id())
+#   # evaluate the density and compare with R
+#   expect_equal(
+#     dmultinom_$evaluate(),
+#     stats::dmultinom(x = x, prob = p, log = TRUE)
+#   )
+#
+#   clear()
+# })
 
-test_that("multinomial_lpdf", {
-  # generate data using R stats:rnorm
-  set.seed(123)
-  p <- (1:10) / sum(1:10)
-  x <- stats::rmultinom(1, 100, p)
-
-  # create a fims Rcpp object
-  # initialize the Dmultinom module
-  dmultinom_ <- new(TMBDmultinomDistribution)
-  # populate class members
-  dmultinom_$x <- new(ParameterVector, x, 10)
-  dmultinom_$expected_values <- new(ParameterVector, p, 10)
-  dmultinom_$is_na <- FALSE
-  dmultinom_$dims <- c(1,10)
-  # evaluate the density and compare with R
-  expect_equal(
-    dmultinom_$evaluate(),
-    stats::dmultinom(x = x, prob = p, log = TRUE)
-  )
-
-  clear()
-})
-
-
-#test something like this:
-
-  test_that("Fleet: SetAgeCompLikelihood works", {
-    fleet <- new(Fleet)
-
-    expect_silent(fleet$SetAgeCompLikelihood(1))
-
-    clear()
-  })
-
-test_that("Fleet: SetIndexLikelihood works", {
-  fleet <- new(Fleet)
 
-  expect_silent(fleet$SetIndexLikelihood(1))
-
-  clear()
-})
-
-test_that("Fleet: SetObservedAgeCompData works", {
-  fleet <- new(Fleet)
-
-  expect_silent(fleet$SetObservedAgeCompData(1))
-
-  clear()
-})
-
-test_that("Fleet: SetObservedIndexData works", {
-  fleet <- new(Fleet)
-
-  expect_silent(fleet$SetObservedIndexData(1))
-
-  clear()
-})