From 3fef7b95613619b990af2ecf3071c18b5eb2a406 Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA <85530309+Andrea-Havron-NOAA@users.noreply.github.com> Date: Fri, 3 Nov 2023 22:47:13 +0000 Subject: [PATCH 1/6] refactor recruitment deviations to log space Co-authored-by: iantaylor-NOAA Co-authored-by: KyleShertzer-NOAA --- inst/include/common/model.hpp | 12 +++---- .../include/interface/rcpp/rcpp_interface.hpp | 6 ++-- .../rcpp/rcpp_objects/rcpp_recruitment.hpp | 32 +++++++++---------- .../population/population.hpp | 8 ++--- .../recruitment/functors/recruitment_base.hpp | 22 ++++++------- tests/gtest/test_population_Recruitment.cpp | 2 +- ...t_population_dynamics_recruitment_base.cpp | 10 +++--- tests/gtest/test_population_test_fixture.hpp | 6 ++-- tests/integration/integration_class.hpp | 8 ++--- ...ration_test_population_tmb_nointerface.cpp | 4 +-- tests/testthat/test-fims-estimation.R | 28 +++++++--------- tests/testthat/test-recruitment-interface.R | 10 +++--- vignettes/fims-demo.Rmd | 22 ++++++------- 13 files changed, 83 insertions(+), 87 deletions(-) diff --git a/inst/include/common/model.hpp b/inst/include/common/model.hpp index e1130014..5d7978ec 100644 --- a/inst/include/common/model.hpp +++ b/inst/include/common/model.hpp @@ -83,7 +83,7 @@ class Model { // may need singleton vector > naa(n_pops); vector > ssb(n_pops); vector > biomass(n_pops); - vector > rec_dev(n_pops); + vector > log_recruit_dev(n_pops); vector > recruitment(n_pops); vector > M(n_pops); #endif @@ -148,8 +148,8 @@ class Model { // may need singleton #ifdef TMB_MODEL naa(pop_idx) = vector((*it).second->numbers_at_age); ssb(pop_idx) = vector((*it).second->spawning_biomass); - rec_dev(pop_idx) = - vector((*it).second->recruitment->recruit_deviations); + log_recruit_dev(pop_idx) = + vector((*it).second->recruitment->log_recruit_devs); recruitment(pop_idx) = vector((*it).second->expected_recruitment); biomass(pop_idx) = vector((*it).second->biomass); M(pop_idx) = vector((*it).second->M); @@ -181,7 +181,7 @@ class Model { // may need singleton REPORT_F(jnll, of); REPORT_F(naa, of); REPORT_F(ssb, of); - REPORT_F(rec_dev, of); + REPORT_F(log_recruit_dev, of); REPORT_F(recruitment, of); REPORT_F(biomass, of); REPORT_F(M, of); @@ -198,7 +198,7 @@ class Model { // may need singleton vector NAA = ADREPORTvector(naa); vector Biomass = ADREPORTvector(biomass); vector SSB = ADREPORTvector(ssb); - vector RecDev = ADREPORTvector(rec_dev); + vector LogRecDev = ADREPORTvector(log_recruit_dev); vector FMort = ADREPORTvector(F_mort); vector ExpectedIndex = ADREPORTvector(exp_index); vector CNAA = ADREPORTvector(cnaa); @@ -206,7 +206,7 @@ class Model { // may need singleton ADREPORT_F(NAA, of); ADREPORT_F(Biomass, of); ADREPORT_F(SSB, of); - ADREPORT_F(RecDev, of); + ADREPORT_F(LogRecDev, of); ADREPORT_F(FMort, of); ADREPORT_F(ExpectedIndex, of); ADREPORT_F(CNAA, of); diff --git a/inst/include/interface/rcpp/rcpp_interface.hpp b/inst/include/interface/rcpp/rcpp_interface.hpp index 45711e9d..ea1d7f91 100644 --- a/inst/include/interface/rcpp/rcpp_interface.hpp +++ b/inst/include/interface/rcpp/rcpp_interface.hpp @@ -350,9 +350,9 @@ RCPP_MODULE(fims) { .constructor() .field("logit_steep", &BevertonHoltRecruitmentInterface::logit_steep) .field("log_rzero", &BevertonHoltRecruitmentInterface::log_rzero) - .field("deviations", &BevertonHoltRecruitmentInterface::deviations) - .field("estimate_deviations", - &BevertonHoltRecruitmentInterface::estimate_deviations) + .field("log_devs", &BevertonHoltRecruitmentInterface::log_devs) + .field("estimate_log_devs", + &BevertonHoltRecruitmentInterface::estimate_log_devs) .method("get_id", &BevertonHoltRecruitmentInterface::get_id) .field("log_sigma_recruit", &BevertonHoltRecruitmentInterface::log_sigma_recruit) diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp index 8eb351d3..713ffeb5 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp @@ -28,7 +28,7 @@ class RecruitmentInterfaceBase : public FIMSRcppInterfaceBase { static std::map live_objects; /**< map associating the ids of RecruitmentInterfaceBase to the objects */ - // static std::vector recruit_deviations; /**< vector of recruitment + // static std::vector log_recruit_devs; /**< vector of log recruitment // deviations*/ /// static bool constrain_deviations; /**< whether or not the rec devs are /// constrained*/ @@ -72,9 +72,9 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { public: Parameter logit_steep; /**< steepness or the productivity of the stock*/ Parameter log_rzero; /**< recruitment at unfished biomass */ - Parameter log_sigma_recruit; /**< the log of the stock recruit deviations */ - Rcpp::NumericVector deviations; /**< recruitment deviations*/ - bool estimate_deviations = + Parameter log_sigma_recruit; /**< the log of the stock recruit standard deviation */ + Rcpp::NumericVector log_devs; /**< log recruitment deviations*/ + bool estimate_log_devs = false; /**< boolean describing whether to estimate */ BevertonHoltRecruitmentInterface() : RecruitmentInterfaceBase() {} @@ -102,13 +102,13 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { fims_popdy::SRBevertonHolt NLL; NLL.log_sigma_recruit = this->log_sigma_recruit.value_m; - NLL.recruit_deviations.resize(deviations.size()); // Vector from TMB - for (int i = 0; i < deviations.size(); i++) { - NLL.recruit_deviations[i] = deviations[i]; + NLL.log_recruit_devs.resize(log_devs.size()); // Vector from TMB + for (int i = 0; i < log_devs.size(); i++) { + NLL.log_recruit_devs[i] = log_devs[i]; } - RECRUITMENT_LOG << "Rec devs being passed to C++ are " << deviations + RECRUITMENT_LOG << "Log recruit devs being passed to C++ are " << log_devs << std::endl; - NLL.estimate_recruit_deviations = this->estimate_deviations; + NLL.estimate_log_recruit_devs = this->estimate_log_devs; return NLL.evaluate_nll(); } @@ -149,15 +149,15 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { } } - recruitment->recruit_deviations.resize(this->deviations.size()); - if (this->estimate_deviations) { - for (size_t i = 0; i < recruitment->recruit_deviations.size(); i++) { - recruitment->recruit_deviations[i] = this->deviations[i]; - info->RegisterParameter(recruitment->recruit_deviations[i]); + recruitment->log_recruit_devs.resize(this->log_devs.size()); + if (this->estimate_log_devs) { + for (size_t i = 0; i < recruitment->log_recruit_devs.size(); i++) { + recruitment->log_recruit_devs[i] = this->log_devs[i]; + info->RegisterParameter(recruitment->log_recruit_devs[i]); } } else { - for (size_t i = 0; i < recruitment->recruit_deviations.size(); i++) { - recruitment->recruit_deviations[i] = this->deviations[i]; + for (size_t i = 0; i < recruitment->log_recruit_devs.size(); i++) { + recruitment->log_recruit_devs[i] = this->log_devs[i]; } } diff --git a/inst/include/population_dynamics/population/population.hpp b/inst/include/population_dynamics/population/population.hpp index 434dccef..da0c35f9 100644 --- a/inst/include/population_dynamics/population/population.hpp +++ b/inst/include/population_dynamics/population/population.hpp @@ -385,8 +385,8 @@ struct Population : public fims_model_object::FIMSObject { POPULATION_LOG << "phi0 = " << phi0 << std::endl; POPULATION_LOG << "spawning_biomass[year - 1] = " << this->spawning_biomass[year - 1] << std::endl; - POPULATION_LOG << "rec devs = " - << this->recruitment->recruit_deviations[year - 1] + POPULATION_LOG << "log recruit devs = " + << this->recruitment->log_recruit_devs[year - 1] << std::endl; POPULATION_LOG << "rec eval = " << this->recruitment->evaluate( @@ -394,9 +394,9 @@ struct Population : public fims_model_object::FIMSObject { << std::endl; this->numbers_at_age[i_age_year] = this->recruitment->evaluate(this->spawning_biomass[year - 1], phi0) * - this->recruitment->recruit_deviations[year]; + fims_math::exp(this->recruitment->log_recruit_devs[year]); this->expected_recruitment[year] = this->numbers_at_age[i_age_year]; - POPULATION_LOG << " numbers at age at indexya " << i_age_year << " is " + POPULATION_LOG << " numbers at age at index i_age_year " << i_age_year << " is " << this->numbers_at_age[i_age_year] << std::endl; } diff --git a/inst/include/population_dynamics/recruitment/functors/recruitment_base.hpp b/inst/include/population_dynamics/recruitment/functors/recruitment_base.hpp index 7e5f096f..62b3ecf1 100644 --- a/inst/include/population_dynamics/recruitment/functors/recruitment_base.hpp +++ b/inst/include/population_dynamics/recruitment/functors/recruitment_base.hpp @@ -32,7 +32,7 @@ struct RecruitmentBase : public fims_model_object::FIMSObject { static uint32_t id_g; /*!< reference id for recruitment object*/ typename fims::ModelTraits::ParameterVector - recruit_deviations; /*!< A vector of recruitment deviations */ + log_recruit_devs; /*!< A vector of recruitment deviations */ bool constrain_deviations = false; /*!< A flag to indicate if recruitment deviations are summing to zero or not */ @@ -40,7 +40,7 @@ struct RecruitmentBase : public fims_model_object::FIMSObject { deviations */ Type log_rzero; /*!< Log of unexploited recruitment.*/ - bool estimate_recruit_deviations = + bool estimate_log_recruit_devs = true; /*!< A flag to indicate if recruitment deviations are estimated or not */ @@ -72,14 +72,14 @@ struct RecruitmentBase : public fims_model_object::FIMSObject { virtual const Type evaluate_nll() { Type nll = 0.0; /*!< The negative log likelihood value */ - if (!this->estimate_recruit_deviations) { + if (!this->estimate_log_recruit_devs) { return nll; } else { #ifdef TMB_MODEL fims_distributions::Dnorm dnorm; dnorm.sd = fims_math::exp(this->log_sigma_recruit); - for (size_t i = 0; i < this->recruit_deviations.size(); i++) { - dnorm.x = fims_math::log(this->recruit_deviations[i]); + for (size_t i = 0; i < this->log_recruit_devs.size(); i++) { + dnorm.x = this->log_recruit_devs[i]; dnorm.mean = 0.0; nll -= dnorm.evaluate(true); } @@ -101,14 +101,14 @@ struct RecruitmentBase : public fims_model_object::FIMSObject { Type sum = 0.0; - for (size_t i = 0; i < this->recruit_deviations.size(); i++) { - sum += this->recruit_deviations[i]; + for (size_t i = 0; i < this->log_recruit_devs.size(); i++) { + sum += this->log_recruit_devs[i]; } - RECRUITMENT_LOG << "recruit_deviations: \n"; - for (size_t i = 0; i < this->recruit_deviations.size(); i++) { - this->recruit_deviations[i] -= sum / (this->recruit_deviations.size()); - RECRUITMENT_LOG << this->recruit_deviations[i] << std::endl; + RECRUITMENT_LOG << "log_recruit_devs: \n"; + for (size_t i = 0; i < this->log_recruit_devs.size(); i++) { + this->log_recruit_devs[i] -= sum / (this->log_recruit_devs.size()); + RECRUITMENT_LOG << this->log_recruit_devs[i] << std::endl; } } }; diff --git a/tests/gtest/test_population_Recruitment.cpp b/tests/gtest/test_population_Recruitment.cpp index d09c71e2..5eaee36e 100644 --- a/tests/gtest/test_population_Recruitment.cpp +++ b/tests/gtest/test_population_Recruitment.cpp @@ -41,7 +41,7 @@ namespace expect_recruitment[r_i_age_year] = (0.8 * rzero * steep * population.spawning_biomass[sb_year]) / - (0.2 * phi0 * rzero * (1.0 - steep) + population.spawning_biomass[sb_year] * (steep - 0.2)) * population.recruitment->recruit_deviations[r_year]; + (0.2 * phi0 * rzero * (1.0 - steep) + population.spawning_biomass[sb_year] * (steep - 0.2)) * fims_math::exp(population.recruitment->log_recruit_devs[r_year]); // calculate recruitment in population module population.CalculateRecruitment(r_i_age_year, r_year); diff --git a/tests/gtest/test_population_dynamics_recruitment_base.cpp b/tests/gtest/test_population_dynamics_recruitment_base.cpp index b37630c5..9961f9bc 100644 --- a/tests/gtest/test_population_dynamics_recruitment_base.cpp +++ b/tests/gtest/test_population_dynamics_recruitment_base.cpp @@ -7,16 +7,16 @@ namespace TEST(RecruitmentDeviations, ConstraintWorks) { fims_popdy::SRBevertonHolt recruit; - recruit.recruit_deviations = {-1.0, 2.0, 3.0}; + recruit.log_recruit_devs = {-1.0, 2.0, 3.0}; // Test if constrain_deviations = false works recruit.constrain_deviations = false; recruit.PrepareConstrainedDeviations(); std::vector expected_deviations_false = {-1.0, 2.0, 3.0}; - for (int i = 0; i < recruit.recruit_deviations.size(); i++) + for (int i = 0; i < recruit.log_recruit_devs.size(); i++) { - EXPECT_EQ(recruit.recruit_deviations[i], + EXPECT_EQ(recruit.log_recruit_devs[i], expected_deviations_false[i]); } @@ -26,9 +26,9 @@ namespace // c(-1.0, 2.0, 3.0)-sum(c(-1.0, 2.0, 3.0))/3 = -2.3333333 0.6666667 1.6666667 std::vector expected_deviations_true = {-2.3333333, 0.6666667, 1.6666667}; - for (int i = 0; i < recruit.recruit_deviations.size(); i++) + for (int i = 0; i < recruit.log_recruit_devs.size(); i++) { - EXPECT_NEAR(recruit.recruit_deviations[i], + EXPECT_NEAR(recruit.log_recruit_devs[i], expected_deviations_true[i], 0.0000001); } } diff --git a/tests/gtest/test_population_test_fixture.hpp b/tests/gtest/test_population_test_fixture.hpp index 0b91e7fc..c8365a6d 100644 --- a/tests/gtest/test_population_test_fixture.hpp +++ b/tests/gtest/test_population_test_fixture.hpp @@ -153,9 +153,9 @@ class PopulationPrepareTestFixture : public testing::Test { auto recruitment = std::make_shared>(); recruitment->logit_steep = fims_math::logit(0.2, 1.0, 0.75); recruitment->log_rzero = fims_math::log(1000000.0); - recruitment->recruit_deviations.resize(nyears); - for (int i = 0; i < recruitment->recruit_deviations.size(); i++) { - recruitment->recruit_deviations[i] = 1.0; + recruitment->log_recruit_devs.resize(nyears); + for (int i = 0; i < recruitment->log_recruit_devs.size(); i++) { + recruitment->log_recruit_devs[i] = 0.0; } population.recruitment = recruitment; } diff --git a/tests/integration/integration_class.hpp b/tests/integration/integration_class.hpp index b89140c7..9ffeb3d7 100644 --- a/tests/integration/integration_class.hpp +++ b/tests/integration/integration_class.hpp @@ -548,8 +548,8 @@ class IntegrationTest { it = obj.find("logR.resid"); - rec->recruit_deviations.resize(nyears + 1); - std::fill(rec->recruit_deviations.begin(), rec->recruit_deviations.end(), 1.0); + rec->log_recruit_devs.resize(nyears + 1); + std::fill(rec->log_recruit_devs.begin(), rec->log_recruit_devs.end(), 0.0); if (it != obj.end()) { if ((*it).second.GetType() == JsonValueType::Array) { JsonArray rdev = (*it).second.GetArray(); @@ -557,9 +557,9 @@ class IntegrationTest { std::cout << "recruitment deviations: "; } for (size_t i = 0; i < rdev.size(); i++) { - rec->recruit_deviations[i] = std::exp(rdev[i].GetDouble()); + rec->log_recruit_devs[i] = rdev[i].GetDouble(); if (print_statements) { - std::cout << rec->recruit_deviations[i] << " "; + std::cout << rec->log_recruit_devs[i] << " "; } } if (print_statements) { diff --git a/tests/integration/integration_test_population_tmb_nointerface.cpp b/tests/integration/integration_test_population_tmb_nointerface.cpp index 9637d5c8..ff5a1e4b 100644 --- a/tests/integration/integration_test_population_tmb_nointerface.cpp +++ b/tests/integration/integration_test_population_tmb_nointerface.cpp @@ -137,8 +137,8 @@ Type objective_function::operator()(){ rec->rzero = R0; rec->steep = h; rec->log_sigma_recruit = logR_sd; - rec->recruit_deviations.resize(nyears); - std::fill(rec->recruit_deviations.begin(), rec->recruit_deviations.end(), 1.0); + rec->log_recruit_devs.resize(nyears); + std::fill(rec->log_recruit_devs.begin(), rec->log_recruit_devs.end(), 0.0); pop.recruitment = rec; //Set maturity diff --git a/tests/testthat/test-fims-estimation.R b/tests/testthat/test-fims-estimation.R index 20344b30..0282ade5 100644 --- a/tests/testthat/test-fims-estimation.R +++ b/tests/testthat/test-fims-estimation.R @@ -38,11 +38,9 @@ setup_fims <- function(om_input, om_output, em_input) { test_env$recruitment$logit_steep$value <- -log(1.0 - om_input$h) + log(om_input$h - 0.2) test_env$recruitment$logit_steep$is_random_effect <- FALSE test_env$recruitment$logit_steep$estimated <- FALSE - test_env$recruitment$estimate_deviations <- TRUE - # recruit deviations should enter the model in normal space. - # The log is taken in the likelihood calculations - # alternative setting: recruitment$deviations <- rep(1, length(om_input$logR.resid)) - test_env$recruitment$deviations <- exp(om_input$logR.resid) + test_env$recruitment$estimate_log_devs <- TRUE + # alternative setting: recruitment$log_devs <- rep(0, length(om_input$logR.resid)) + test_env$recruitment$log_devs <- om_input$logR.resid # Data test_env$catch <- em_input$L.obs$fleet1 @@ -209,8 +207,8 @@ test_that("deterministic test of fims", { expect_equal(report$recruitment[[1]][i], om_output$N.age[i, 1]) } - # recruitment deviations (fixed at initial "true" values) - expect_equal(log(report$rec_dev[[1]]), om_input$logR.resid) + # recruitment log_devs (fixed at initial "true" values) + expect_equal(report$log_recruit_dev[[1]], om_input$logR.resid) # F (fixed at initial "true" values) expect_equal(report$F_mort[[1]], om_output$f) @@ -312,7 +310,7 @@ test_that("nll test of fims", { # recruitment likelihood rec_nll <- -sum(dnorm( - log(nll_env$recruitment$deviations), rep(0, om_input$nyr), + nll_env$recruitment$log_devs, rep(0, om_input$nyr), om_input$logR_sd, TRUE )) @@ -440,12 +438,12 @@ test_that("estimation test of fims", { report$recruitment[[1]][1:om_input$nyr] ) - # recruitment deviations - sdr_rdev <- sdr_report[which(rownames(sdr_report) == "RecDev"), ] + # recruitment log deviations + sdr_rdev <- sdr_report[which(rownames(sdr_report) == "LogRecDev"), ] rdev_are <- rep(0, length(om_input$logR.resid)) for (i in 1:length(om_input$logR.resid)) { - rdev_are[i] <- abs(report$rec_dev[[1]][i] - exp(om_input$logR.resid[i])) # / + rdev_are[i] <- abs(report$log_recruit_dev[[1]][i] - om_input$logR.resid[i]) # / # exp(om_input$logR.resid[i]) # expect_lte(rdev_are[i], 1) # 1 } @@ -585,11 +583,9 @@ test_that("run FIMS in a for loop", { recruitment$logit_steep$value <- -log(1.0 - om_input$h) + log(om_input$h - 0.2) recruitment$logit_steep$is_random_effect <- FALSE recruitment$logit_steep$estimated <- FALSE - recruitment$estimate_deviations <- TRUE - # recruit deviations should enter the model in normal space. - # The log is taken in the likelihood calculations - # alternative setting: recruitment$deviations <- rep(1, length(om_input$logR.resid)) - recruitment$deviations <- exp(om_input$logR.resid) + recruitment$estimate_log_devs <- TRUE + # alternative setting: recruitment$log_devs <- rep(0, length(om_input$logR.resid)) + recruitment$log_devs <- om_input$logR.resid # Data catch <- em_input$L.obs$fleet1 diff --git a/tests/testthat/test-recruitment-interface.R b/tests/testthat/test-recruitment-interface.R index e95e5bc8..22f41e83 100644 --- a/tests/testthat/test-recruitment-interface.R +++ b/tests/testthat/test-recruitment-interface.R @@ -29,17 +29,17 @@ test_that("Recruitment input settings work as expected", { expect_equal(object = recruitment$evaluate(spawns, ssb0), expected = 1090802.68) - devs <- c(1.0, 2.0, 3.0) - recruitment$deviations <- devs + log_devs <- c(-1.0, 2.0, 3.0) + recruitment$log_devs <- log_devs - expected_nll <- -sum(log(stats::dnorm(log(devs), 0, 0.7))) + expected_nll <- -sum(log(stats::dnorm(log_devs, 0, 0.7))) - recruitment$estimate_deviations <- FALSE + recruitment$estimate_log_devs <- FALSE expect_equal(recruitment$evaluate_nll(), 0.0) - recruitment$estimate_deviations <- TRUE + recruitment$estimate_log_devs <- TRUE expect_equal(recruitment$evaluate_nll(), expected = expected_nll) fims$clear() diff --git a/vignettes/fims-demo.Rmd b/vignettes/fims-demo.Rmd index 0c093cee..58658957 100644 --- a/vignettes/fims-demo.Rmd +++ b/vignettes/fims-demo.Rmd @@ -252,19 +252,19 @@ recruitment$logit_steep$is_random_effect <- FALSE recruitment$logit_steep$estimated <- FALSE ``` -We also need to set up recruitment deviations. FIMS recruitment modules have a boolean, *estimate_deviations* to specify whether or not deviations are estimated; and a vector, *deviations* to set the deviation values. The *deviations* vector takes values greater than 0 as these will be logged in the likelihood. +We also need to set up log recruitment deviations. FIMS recruitment modules have a boolean, *estimate_log_devs* to specify whether or not log deviations are estimated; and a vector, *log_devs* to set the log deviation values. ```{r rec-devs} -recruitment$estimate_deviations <- FALSE -recruitment$deviations <- c( - 1.0931337, 1.5494153, 0.8754735, 0.6488721, - 1.9128839, 1.6593211, 0.9327825, 1.3531871, - 0.9207434, 1.2304792, 1.1652038, 0.8048559, - 0.8752845, 1.1187967, 0.8989675, 1.3083559, - 1.2724463, 0.5799534, 0.7891447, 0.5571984, - 1.3515173, 1.2452116, 0.8002613, 0.5983474, - 1.1704664, 0.5828168, 0.8223697, 1.2225558, - 1.4513402, 0.9308739 +recruitment$estimate_log_devs <- FALSE +recruitment$log_devs <- c( + 0.08904850, 0.43787763, -0.13299042, -0.43251973, + 0.64861200, 0.50640852, -0.06958319, 0.30246260, + -0.08257384, 0.20740372, 0.15289604, -0.21709207, + -0.13320626, 0.11225374, -0.10650836, 0.26877132, + 0.24094126, -0.54480751, -0.23680557, -0.58483386, + 0.30122785, 0.21930545, -0.22281699, -0.51358369, + 0.15740234, -0.53988240, -0.19556523, 0.20094360, + 0.37248740, -0.07163145 ) ``` #### Growth From 67cd1a97f9bda8d27ac1baecdedf85b97b494f09 Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA <85530309+Andrea-Havron-NOAA@users.noreply.github.com> Date: Mon, 6 Nov 2023 19:58:31 +0000 Subject: [PATCH 2/6] fix rec dev vector length --- .../population_dynamics/population/population.hpp | 2 +- tests/testthat/test-fims-estimation.R | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/inst/include/population_dynamics/population/population.hpp b/inst/include/population_dynamics/population/population.hpp index da0c35f9..4710674e 100644 --- a/inst/include/population_dynamics/population/population.hpp +++ b/inst/include/population_dynamics/population/population.hpp @@ -394,7 +394,7 @@ struct Population : public fims_model_object::FIMSObject { << std::endl; this->numbers_at_age[i_age_year] = this->recruitment->evaluate(this->spawning_biomass[year - 1], phi0) * - fims_math::exp(this->recruitment->log_recruit_devs[year]); + fims_math::exp(this->recruitment->log_recruit_devs[year-1]); this->expected_recruitment[year] = this->numbers_at_age[i_age_year]; POPULATION_LOG << " numbers at age at index i_age_year " << i_age_year << " is " << this->numbers_at_age[i_age_year] << std::endl; diff --git a/tests/testthat/test-fims-estimation.R b/tests/testthat/test-fims-estimation.R index 0282ade5..0aa7ce74 100644 --- a/tests/testthat/test-fims-estimation.R +++ b/tests/testthat/test-fims-estimation.R @@ -577,15 +577,15 @@ test_that("run FIMS in a for loop", { # logR_sd is NOT logged. It needs to enter the model logged b/c the exp() is taken # before the likelihood calculation recruitment$log_sigma_recruit$value <- log(om_input$logR_sd) - recruitment$log_rzero$value <- log(om_input$R0) + recruitment$log_rzero$value <- 13 #log(om_input$R0) recruitment$log_rzero$is_random_effect <- FALSE recruitment$log_rzero$estimated <- TRUE recruitment$logit_steep$value <- -log(1.0 - om_input$h) + log(om_input$h - 0.2) recruitment$logit_steep$is_random_effect <- FALSE recruitment$logit_steep$estimated <- FALSE recruitment$estimate_log_devs <- TRUE - # alternative setting: recruitment$log_devs <- rep(0, length(om_input$logR.resid)) - recruitment$log_devs <- om_input$logR.resid + recruitment$log_devs <- rep(0, length(om_input$logR.resid)-1) + #recruitment$log_devs <- om_input$logR.resid[-1] # Data catch <- em_input$L.obs$fleet1 @@ -692,16 +692,16 @@ test_that("run FIMS in a for loop", { fims$CreateTMBModel() parameters <- list(p = 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) )) - + report <- obj$report(obj$par) expect_false(is.null(report)) - max_gradient <- max(abs(obj$gr(opt$par))) + max_gradient <- max(abs(obj$gr(obj$env$last.par.best))) expect_lte(max_gradient, 0.00001) fims$clear() } From 9f1ff77939e06a106feb0598b96aad40368c6145 Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA <85530309+Andrea-Havron-NOAA@users.noreply.github.com> Date: Wed, 8 Nov 2023 01:03:28 +0000 Subject: [PATCH 3/6] fix C++ tests --- .../population_dynamics/population/population.hpp | 12 ++++++++---- tests/gtest/test_population_Recruitment.cpp | 7 +++++-- tests/integration/integration_class.hpp | 9 ++++++--- .../integration_test_population_tmb_nointerface.cpp | 5 ++++- 4 files changed, 23 insertions(+), 10 deletions(-) diff --git a/inst/include/population_dynamics/population/population.hpp b/inst/include/population_dynamics/population/population.hpp index 4710674e..684a5b0b 100644 --- a/inst/include/population_dynamics/population/population.hpp +++ b/inst/include/population_dynamics/population/population.hpp @@ -377,8 +377,9 @@ struct Population : public fims_model_object::FIMSObject { * * @param i_age_year dimension folded index for age and year * @param year the year recruitment is being calculated for + * @param i_dev index to log_recruit_dev of vector length nyears-1 */ - void CalculateRecruitment(size_t i_age_year, size_t year) { + void CalculateRecruitment(size_t i_age_year, size_t year, size_t i_dev) { POPULATION_LOG << "recruitment 1" << std::endl; Type phi0 = CalculateSBPR0(); POPULATION_LOG << "recruitment 2" << std::endl; @@ -386,7 +387,7 @@ struct Population : public fims_model_object::FIMSObject { POPULATION_LOG << "spawning_biomass[year - 1] = " << this->spawning_biomass[year - 1] << std::endl; POPULATION_LOG << "log recruit devs = " - << this->recruitment->log_recruit_devs[year - 1] + << this->recruitment->log_recruit_devs[i_dev - 1] << std::endl; POPULATION_LOG << "rec eval = " << this->recruitment->evaluate( @@ -394,7 +395,10 @@ struct Population : public fims_model_object::FIMSObject { << std::endl; this->numbers_at_age[i_age_year] = this->recruitment->evaluate(this->spawning_biomass[year - 1], phi0) * - fims_math::exp(this->recruitment->log_recruit_devs[year-1]); + /*the log_recruit_dev vector does not include a value for year == 0 + and is of length nyears - 1 where the first position of the vector + corresponds to the second year of the time series.*/ + fims_math::exp(this->recruitment->log_recruit_devs[i_dev-1]); this->expected_recruitment[year] = this->numbers_at_age[i_age_year]; POPULATION_LOG << " numbers at age at index i_age_year " << i_age_year << " is " << this->numbers_at_age[i_age_year] << std::endl; @@ -644,7 +648,7 @@ struct Population : public fims_model_object::FIMSObject { // Set the nrecruits for age a=0 year y (use pointers instead of // functional returns) assuming fecundity = 1 and 50:50 sex ratio POPULATION_LOG << "Recruitment: " << std::endl; - CalculateRecruitment(i_age_year, y); + CalculateRecruitment(i_age_year, y, y); this->unfished_numbers_at_age[i_age_year] = fims_math::exp(this->recruitment->log_rzero); diff --git a/tests/gtest/test_population_Recruitment.cpp b/tests/gtest/test_population_Recruitment.cpp index af212efa..80c89e13 100644 --- a/tests/gtest/test_population_Recruitment.cpp +++ b/tests/gtest/test_population_Recruitment.cpp @@ -41,10 +41,13 @@ namespace expect_recruitment[r_i_age_year] = (0.8 * rzero * steep * population.spawning_biomass[sb_year]) / - (0.2 * phi0 * rzero * (1.0 - steep) + population.spawning_biomass[sb_year] * (steep - 0.2)) * fims_math::exp(population.recruitment->log_recruit_devs[r_year]); + /*the log_recruit_dev vector does not include a value for year == 0 + and is of length nyears - 1 where the first position of the vector + corresponds to the second year of the time series.*/ + (0.2 * phi0 * rzero * (1.0 - steep) + population.spawning_biomass[sb_year] * (steep - 0.2)) * fims_math::exp(population.recruitment->log_recruit_devs[r_year-1]); // calculate recruitment in population module - population.CalculateRecruitment(r_i_age_year, r_year); + population.CalculateRecruitment(r_i_age_year, r_year, r_year); // testing that expected recruitment and population.numbers_at_age match // EXPECT_DOUBLE_EQ() verifies that the two double values are approximately equal, to within 4 ULPs from each other. diff --git a/tests/integration/integration_class.hpp b/tests/integration/integration_class.hpp index 9ffeb3d7..543330a6 100644 --- a/tests/integration/integration_class.hpp +++ b/tests/integration/integration_class.hpp @@ -548,7 +548,10 @@ class IntegrationTest { it = obj.find("logR.resid"); - rec->log_recruit_devs.resize(nyears + 1); + /*the log_recruit_dev vector does not include a value for year == 0 + and is of length nyears - 1 where the first position of the vector + corresponds to the second year of the time series.*/ + rec->log_recruit_devs.resize(nyears); std::fill(rec->log_recruit_devs.begin(), rec->log_recruit_devs.end(), 0.0); if (it != obj.end()) { if ((*it).second.GetType() == JsonValueType::Array) { @@ -556,8 +559,8 @@ class IntegrationTest { if (print_statements) { std::cout << "recruitment deviations: "; } - for (size_t i = 0; i < rdev.size(); i++) { - rec->log_recruit_devs[i] = rdev[i].GetDouble(); + for (size_t i = 0; i < rec->log_recruit_devs.size(); i++) { + rec->log_recruit_devs[i] = rdev[i+1].GetDouble(); if (print_statements) { std::cout << rec->log_recruit_devs[i] << " "; } diff --git a/tests/integration/integration_test_population_tmb_nointerface.cpp b/tests/integration/integration_test_population_tmb_nointerface.cpp index ff5a1e4b..46b19b66 100644 --- a/tests/integration/integration_test_population_tmb_nointerface.cpp +++ b/tests/integration/integration_test_population_tmb_nointerface.cpp @@ -137,7 +137,10 @@ Type objective_function::operator()(){ rec->rzero = R0; rec->steep = h; rec->log_sigma_recruit = logR_sd; - rec->log_recruit_devs.resize(nyears); + /*the log_recruit_dev vector does not include a value for year == 0 + and is of length nyears - 1 where the first position of the vector + corresponds to the second year of the time series.*/ + rec->log_recruit_devs.resize(nyears-1); std::fill(rec->log_recruit_devs.begin(), rec->log_recruit_devs.end(), 0.0); pop.recruitment = rec; From 6c01f126e5d0ad1d7bec23fbf5099383f620098b Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA <85530309+Andrea-Havron-NOAA@users.noreply.github.com> Date: Wed, 8 Nov 2023 01:41:39 +0000 Subject: [PATCH 4/6] fix R tests --- tests/testthat/test-fims-estimation.R | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/tests/testthat/test-fims-estimation.R b/tests/testthat/test-fims-estimation.R index 0aa7ce74..7164edf1 100644 --- a/tests/testthat/test-fims-estimation.R +++ b/tests/testthat/test-fims-estimation.R @@ -40,7 +40,7 @@ setup_fims <- function(om_input, om_output, em_input) { test_env$recruitment$logit_steep$estimated <- FALSE test_env$recruitment$estimate_log_devs <- TRUE # alternative setting: recruitment$log_devs <- rep(0, length(om_input$logR.resid)) - test_env$recruitment$log_devs <- om_input$logR.resid + test_env$recruitment$log_devs <- om_input$logR.resid[-1] # Data test_env$catch <- em_input$L.obs$fleet1 @@ -208,7 +208,8 @@ test_that("deterministic test of fims", { } # recruitment log_devs (fixed at initial "true" values) - expect_equal(report$log_recruit_dev[[1]], om_input$logR.resid) + # 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) @@ -309,8 +310,9 @@ test_that("nll test of fims", { 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), + nll_env$recruitment$log_devs, rep(0, om_input$nyr-1), om_input$logR_sd, TRUE )) @@ -439,16 +441,17 @@ test_that("estimation test of fims", { ) # 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)) + rdev_are <- rep(0, length(om_input$logR.resid)-1) - for (i in 1:length(om_input$logR.resid)) { - rdev_are[i] <- abs(report$log_recruit_dev[[1]][i] - om_input$logR.resid[i]) # / + for (i in 1:length(report$log_recruit_dev[[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), 2]), + sum(rdev_are > qnorm(.975) * sdr_rdev[1:length(om_input$logR.resid)-1, 2]), 0.05 * length(om_input$logR.resid) ) @@ -585,7 +588,6 @@ test_that("run FIMS in a for loop", { recruitment$logit_steep$estimated <- FALSE recruitment$estimate_log_devs <- TRUE recruitment$log_devs <- rep(0, length(om_input$logR.resid)-1) - #recruitment$log_devs <- om_input$logR.resid[-1] # Data catch <- em_input$L.obs$fleet1 From 7d4ed6fd6e7a04575a7491dfdf193735b36c5f18 Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA <85530309+Andrea-Havron-NOAA@users.noreply.github.com> Date: Fri, 10 Nov 2023 00:06:35 +0000 Subject: [PATCH 5/6] fix tests --- tests/gtest/test_population_Recruitment.cpp | 4 ++-- tests/gtest/test_population_test_fixture.hpp | 9 ++++++--- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/tests/gtest/test_population_Recruitment.cpp b/tests/gtest/test_population_Recruitment.cpp index 80c89e13..b0fc8645 100644 --- a/tests/gtest/test_population_Recruitment.cpp +++ b/tests/gtest/test_population_Recruitment.cpp @@ -39,11 +39,11 @@ namespace // vector for storing expected recruitment std::vector expect_recruitment(population.nyears * population.nages, 0.0); - expect_recruitment[r_i_age_year] = - (0.8 * rzero * steep * population.spawning_biomass[sb_year]) / /*the log_recruit_dev vector does not include a value for year == 0 and is of length nyears - 1 where the first position of the vector corresponds to the second year of the time series.*/ + expect_recruitment[r_i_age_year] = + (0.8 * rzero * steep * population.spawning_biomass[sb_year]) / (0.2 * phi0 * rzero * (1.0 - steep) + population.spawning_biomass[sb_year] * (steep - 0.2)) * fims_math::exp(population.recruitment->log_recruit_devs[r_year-1]); // calculate recruitment in population module diff --git a/tests/gtest/test_population_test_fixture.hpp b/tests/gtest/test_population_test_fixture.hpp index 283f13d6..e405a950 100644 --- a/tests/gtest/test_population_test_fixture.hpp +++ b/tests/gtest/test_population_test_fixture.hpp @@ -167,9 +167,12 @@ namespace auto recruitment = std::make_shared>(); recruitment->logit_steep = fims_math::logit(0.2, 1.0, 0.75); recruitment->log_rzero = fims_math::log(1000000.0); - recruitment->recruit_deviations.resize(nyears); - for (int i = 0; i < recruitment->recruit_deviations.size(); i++) { - recruitment->recruit_deviations[i] = 1.0; + /*the log_recruit_dev vector does not include a value for year == 0 + and is of length nyears - 1 where the first position of the vector + corresponds to the second year of the time series.*/ + recruitment->log_recruit_devs.resize(nyears-1); + for (int i = 0; i < recruitment->log_recruit_devs.size(); i++) { + recruitment->log_recruit_devs[i] = 0.0; } population.recruitment = recruitment; From 4c9105c01a36ce48f02660e41d110873ac37be7e Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA <85530309+Andrea-Havron-NOAA@users.noreply.github.com> Date: Fri, 10 Nov 2023 00:19:08 +0000 Subject: [PATCH 6/6] change tolerance in gradient tests to 1e-4 --- tests/testthat/test-fims-estimation.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-fims-estimation.R b/tests/testthat/test-fims-estimation.R index 7164edf1..1d5a6cdc 100644 --- a/tests/testthat/test-fims-estimation.R +++ b/tests/testthat/test-fims-estimation.R @@ -704,7 +704,7 @@ test_that("run FIMS in a for loop", { expect_false(is.null(report)) max_gradient <- max(abs(obj$gr(obj$env$last.par.best))) - expect_lte(max_gradient, 0.00001) + expect_lte(max_gradient, 0.0001) fims$clear() } })