From 32f89db15631572cc342e58034c1dd7cdc647331 Mon Sep 17 00:00:00 2001 From: Bai-Li-NOAA Date: Wed, 16 Aug 2023 13:45:37 -0400 Subject: [PATCH] Co-authored-by: Chris Legault --- .../functors/tmb_distributions.hpp | 12 ++--- .../include/interface/rcpp/rcpp_interface.hpp | 8 +--- .../rcpp/rcpp_objects/rcpp_recruitment.hpp | 16 ------- .../rcpp_objects/rcpp_tmb_distribution.hpp | 3 -- .../population/population.hpp | 2 +- .../recruitment/functors/recruitment_base.hpp | 45 +---------------- ...t_population_dynamics_recruitment_base.cpp | 48 ------------------- tests/integration/integration_class.hpp | 1 - tests/testthat/test-recruitment-interface.R | 12 ----- tests/testthat/test-tmb-distributions.R | 1 - 10 files changed, 8 insertions(+), 140 deletions(-) diff --git a/inst/include/distributions/functors/tmb_distributions.hpp b/inst/include/distributions/functors/tmb_distributions.hpp index 44a9cc5a6..071a47125 100644 --- a/inst/include/distributions/functors/tmb_distributions.hpp +++ b/inst/include/distributions/functors/tmb_distributions.hpp @@ -77,9 +77,7 @@ struct Dlnorm : public DistributionsBase { T x; /*!< observation */ T meanlog; /*!< mean of the distribution of log(x) */ T sdlog; /*!< standard deviation of the distribution of log(x) */ - bool do_bias_correction = - false; /*!< whether or not to bias correct the distribution */ - + Dlnorm() : DistributionsBase() {} /** @@ -94,11 +92,9 @@ struct Dlnorm : public DistributionsBase { virtual const T evaluate(const bool& do_log) { T logx = log(x); T nll; - if (do_bias_correction) { - nll = dnorm(logx, meanlog - pow(sdlog, 2) / 2, sdlog, true) - logx; - } else { - nll = dnorm(logx, meanlog, sdlog, true) - logx; - } + + nll = dnorm(logx, meanlog, sdlog, true) - logx; + if (do_log) { return nll; } else { diff --git a/inst/include/interface/rcpp/rcpp_interface.hpp b/inst/include/interface/rcpp/rcpp_interface.hpp index ecaf74aaf..7b85c2f3e 100644 --- a/inst/include/interface/rcpp/rcpp_interface.hpp +++ b/inst/include/interface/rcpp/rcpp_interface.hpp @@ -188,10 +188,6 @@ RCPP_MODULE(fims) { .field("estimate_deviations", &BevertonHoltRecruitmentInterface::estimate_deviations) .method("get_id", &BevertonHoltRecruitmentInterface::get_id) - .field("recruitment_bias_adj", - &BevertonHoltRecruitmentInterface::recruit_bias_adjustment) - .field("use_bias_correction", - &BevertonHoltRecruitmentInterface::use_bias_correction) .field("log_sigma_recruit", &BevertonHoltRecruitmentInterface::log_sigma_recruit) .method("evaluate", &BevertonHoltRecruitmentInterface::evaluate) @@ -293,9 +289,7 @@ RCPP_MODULE(fims) { .method("evaluate", &DlnormDistributionsInterface::evaluate) .field("x", &DlnormDistributionsInterface::x) .field("meanlog", &DlnormDistributionsInterface::meanlog) - .field("sdlog", &DlnormDistributionsInterface::sdlog) - .field("do_bias_correction", - &DlnormDistributionsInterface::do_bias_correction); + .field("sdlog", &DlnormDistributionsInterface::sdlog); Rcpp::class_("TMBDmultinomDistribution") .constructor() diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp index c8b6062cd..97b265f7f 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp @@ -31,8 +31,6 @@ class RecruitmentInterfaceBase : public FIMSRcppInterfaceBase { // deviations*/ /// static bool constrain_deviations; /**< whether or not the rec devs are /// constrained*/ - // static std::vector rec_bias_adj; /**< a vector of bias adjustment - // values*/ RecruitmentInterfaceBase() { this->id = RecruitmentInterfaceBase::id_g++; @@ -72,12 +70,9 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { 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 recruit_bias_adjustment; /**log_sigma_recruit.value; NLL.recruit_deviations.resize(deviations.size()); // Vector from TMB - NLL.recruit_bias_adjustment.resize( - recruit_bias_adjustment.size()); // Vector from TMB for (int i = 0; i < deviations.size(); i++) { NLL.recruit_deviations[i] = deviations[i]; - NLL.recruit_bias_adjustment[i] = recruit_bias_adjustment[i]; } FIMS_LOG << "Rec devs being passed to C++ are " << deviations << std::endl; - Rcout << "Rec bias adj being passed to C++ are " << recruit_bias_adjustment - << std::endl; - - NLL.use_recruit_bias_adjustment = this->use_bias_correction; NLL.estimate_recruit_deviations = this->estimate_deviations; return NLL.evaluate_nll(); } @@ -170,7 +158,6 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { } } - b0->use_recruit_bias_adjustment = this->use_bias_correction; // add to Information d0->recruitment_models[b0->id] = b0; @@ -220,7 +207,6 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { } } - b1->use_recruit_bias_adjustment = this->use_bias_correction; // add to Information d1->recruitment_models[b1->id] = b1; @@ -270,7 +256,6 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { } } - b2->use_recruit_bias_adjustment = this->use_bias_correction; // add to Information d2->recruitment_models[b2->id] = b2; @@ -320,7 +305,6 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { } } - b3->use_recruit_bias_adjustment = this->use_bias_correction; // add to Information d3->recruitment_models[b3->id] = b3; diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_tmb_distribution.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_tmb_distribution.hpp index e169f148f..596c3479d 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_tmb_distribution.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_tmb_distribution.hpp @@ -161,8 +161,6 @@ class DlnormDistributionsInterface : public DistributionsInterfaceBase { Parameter x; /*!< observation */ Parameter meanlog; /*!< mean of the distribution of log(x) */ Parameter sdlog; /*!< standard deviation of the distribution of log(x) */ - bool do_bias_correction; /*!< true if the lognormal should be bias corrected, - default FALSE */ DlnormDistributionsInterface() : DistributionsInterfaceBase() {} @@ -185,7 +183,6 @@ class DlnormDistributionsInterface : public DistributionsInterfaceBase { dlnorm.x = this->x.value; dlnorm.meanlog = this->meanlog.value; dlnorm.sdlog = this->sdlog.value; - dlnorm.do_bias_correction = this->do_bias_correction; return dlnorm.evaluate(do_log); } diff --git a/inst/include/population_dynamics/population/population.hpp b/inst/include/population_dynamics/population/population.hpp index b4b8d795b..ca7496797 100644 --- a/inst/include/population_dynamics/population/population.hpp +++ b/inst/include/population_dynamics/population/population.hpp @@ -535,7 +535,7 @@ struct Population : public FIMSObject { /* Sets derived vectors to zero Performs parameters transformations - Sets recruitment deviations to mean 0, and then adds bias adjustment. + Sets recruitment deviations to mean 0. */ Prepare(); /* diff --git a/inst/include/population_dynamics/recruitment/functors/recruitment_base.hpp b/inst/include/population_dynamics/recruitment/functors/recruitment_base.hpp index 5fe24f2be..315875428 100644 --- a/inst/include/population_dynamics/recruitment/functors/recruitment_base.hpp +++ b/inst/include/population_dynamics/recruitment/functors/recruitment_base.hpp @@ -35,18 +35,7 @@ struct RecruitmentBase : public FIMSObject { recruit_deviations; /*!< A vector of recruitment deviations */ bool constrain_deviations = false; /*!< A flag to indicate if recruitment deviations are summing to zero or not */ - typename ModelTraits::DataVector - recruit_bias_adjustment; /*!< A vector of bias adj values - (incorporating sigma_recruit)*/ - // Initially fixing bias adjustment (b_y in collobarative - // workflow specification) to 1.0. - // In the future, this would be set by the user. - typename ModelTraits::DataVector - recruit_bias_adjustment_fraction; /*!< A vector of bias adjustment - fractions (on the 0 to 1 range)*/ - bool use_recruit_bias_adjustment = - true; /*!< A flag to indicate if recruitment deviations are bias adjusted - */ + Type log_sigma_recruit; /*!< Log standard deviation of log recruitment deviations */ Type log_rzero; /*!< Log of unexploited recruitment.*/ @@ -62,13 +51,10 @@ struct RecruitmentBase : public FIMSObject { virtual ~RecruitmentBase() {} /** - * @brief Prepares the recruitment bias adjustment vector. + * @brief Prepares the recruitment deviations vector. * */ void Prepare() { - this->recruit_bias_adjustment_fraction.resize( - this->recruit_deviations.size()); - this->recruit_bias_adjustment.resize(this->recruit_deviations.size()); this->PrepareConstrainedDeviations(); } @@ -97,9 +83,6 @@ struct RecruitmentBase : public FIMSObject { for (size_t i = 0; i < this->recruit_deviations.size(); i++) { dnorm.x = fims::log(this->recruit_deviations[i]); dnorm.mean = 0.0; - if (this->use_recruit_bias_adjustment) { - dnorm.mean -= this->recruit_bias_adjustment[i]; - } nll -= dnorm.evaluate(true); } #endif @@ -131,30 +114,6 @@ struct RecruitmentBase : public FIMSObject { } } - /** @brief Prepare recruitment bias adjustment. - * Based on Methot & Taylor (2011). - */ - void PrepareBiasAdjustment() { - Type recruit_bias_adjustment_size = this->recruit_bias_adjustment.size(); - - if (!this->use_recruit_bias_adjustment) { - for (size_t i = 0; i < recruit_bias_adjustment_size; i++) { - this->recruit_bias_adjustment_fraction[i] = 1.0; - this->recruit_bias_adjustment[i] = 0.0; - } - } else { - for (size_t i = 0; i < recruit_bias_adjustment_size; i++) { - // Initially fixing bias adjustment (b_y in collobarative - // workflow specification) to 1.0. - // In the future, this would be set by the user. - this->recruit_bias_adjustment_fraction[i] = 1.0; - this->recruit_bias_adjustment[i] = - 0.5 * fims::exp(this->log_sigma_recruit) * - fims::exp(this->log_sigma_recruit) * - this->recruit_bias_adjustment_fraction[i]; - } - } - } }; template diff --git a/tests/gtest/test_population_dynamics_recruitment_base.cpp b/tests/gtest/test_population_dynamics_recruitment_base.cpp index 37eecfcda..926601966 100644 --- a/tests/gtest/test_population_dynamics_recruitment_base.cpp +++ b/tests/gtest/test_population_dynamics_recruitment_base.cpp @@ -33,52 +33,4 @@ namespace } } - TEST(recruitment_bias_adjustment, PrepareBiasAdjustment_works) - { - fims::SRBevertonHolt recruit; - recruit.log_sigma_recruit = fims::log(0.7); - - // Test if use_recruit_bias_adjustment = false works - recruit.use_recruit_bias_adjustment = false; - recruit.recruit_bias_adjustment = {0.0, 0.0, 0.0}; - recruit.recruit_bias_adjustment_fraction = {1.0, 1.0, 1.0}; - recruit.PrepareBiasAdjustment(); - - std::vector expected_bias_adjustment_false = {0.0, 0.0, 0.0}; - for (int i = 0; i < recruit.recruit_bias_adjustment.size(); i++) - { - EXPECT_EQ(recruit.recruit_bias_adjustment[i], - expected_bias_adjustment_false[i]); - } - - // Test if use_recruit_bias_adjustment = true works - recruit.use_recruit_bias_adjustment = true; - recruit.recruit_bias_adjustment = {0.0, 0.0, 0.0}; - recruit.recruit_bias_adjustment_fraction = {1.0, 1.0, 1.0}; - recruit.PrepareBiasAdjustment(); - - // R code to generate true values: 0.5 * 0.7^2 * c(1.0, 1.0, 1.0) - // 0.245 0.245 0.245 - std::vector expected_bias_adjustment_true = {0.245, 0.245, 0.245}; - for (int i = 0; i < recruit.recruit_bias_adjustment.size(); i++) - { - EXPECT_NEAR(recruit.recruit_bias_adjustment[i], - expected_bias_adjustment_true[i], 0.001); - } - - // Test if recruit_bias_adjustment_fraction will be converted to 1.0 - // when initial input value is 2.0 - recruit.recruit_bias_adjustment_fraction = {2.0, 2.0, 2.0}; - recruit.PrepareBiasAdjustment(); - - // R code to generate true values: 0.5 * 0.7^2 * c(1.0, 1.0, 1.0) - // 0.245 0.245 0.245 - std::vector expected_bias_adjustment = {0.245, 0.245, 0.245}; - for (int i = 0; i < recruit.recruit_bias_adjustment.size(); i++) - { - EXPECT_NEAR(recruit.recruit_bias_adjustment[i], - expected_bias_adjustment[i], 0.001); - } - } - } diff --git a/tests/integration/integration_class.hpp b/tests/integration/integration_class.hpp index d41fccb5c..947030174 100644 --- a/tests/integration/integration_class.hpp +++ b/tests/integration/integration_class.hpp @@ -571,7 +571,6 @@ class IntegrationTest { std::cout << "'logR.resid' not found.\n"; } } - rec->use_recruit_bias_adjustment = false; pop.recruitment = rec; // set maturity diff --git a/tests/testthat/test-recruitment-interface.R b/tests/testthat/test-recruitment-interface.R index 872ffbcf5..e95e5bc85 100644 --- a/tests/testthat/test-recruitment-interface.R +++ b/tests/testthat/test-recruitment-interface.R @@ -17,8 +17,6 @@ test_that("Recruitment input settings work as expected", { recruitment$logit_steep$estimated <- TRUE recruitment$log_rzero$value <- log(r0) recruitment$log_sigma_recruit$value <- log(0.7) - recruitment$recruitment_bias_adj <- rep(1.0, 3) - recruitment$use_bias_correction <- FALSE expect_equal(recruitment$get_id(), 1) expect_equal(recruitment$logit_steep$value, 0.78845736) @@ -44,15 +42,5 @@ test_that("Recruitment input settings work as expected", { recruitment$estimate_deviations <- TRUE expect_equal(recruitment$evaluate_nll(), expected = expected_nll) - recruitment$use_bias_correction <- TRUE - recruitment$recruitment_bias_adj <- rep(0.245, 3) - expected_nll <- -sum(log(stats::dnorm(log(devs), -0.245, 0.7))) - expect_equal(recruitment$evaluate_nll(), expected = expected_nll) - - recruitment$use_bias_correction <- TRUE - recruitment$recruitment_bias_adj <- c(0.245, 0.2, 0.1) - expected_nll <- -sum(log(stats::dnorm(log(devs), c(-0.245, -0.2, -0.1), 0.7))) - expect_equal(recruitment$evaluate_nll(), expected = expected_nll) - fims$clear() }) diff --git a/tests/testthat/test-tmb-distributions.R b/tests/testthat/test-tmb-distributions.R index 8f84b0ca1..892a2ac0f 100644 --- a/tests/testthat/test-tmb-distributions.R +++ b/tests/testthat/test-tmb-distributions.R @@ -34,7 +34,6 @@ test_that("dlnorm", { dlnorm_$x$value <- y dlnorm_$meanlog$value <- 0 dlnorm_$sdlog$value <- 1 - dlnorm_$do_bias_correction <- FALSE # evaluate the density and compare with R expect_equal(dlnorm_$evaluate(TRUE), stats::dlnorm(y, 0, 1, TRUE)) expect_equal(dlnorm_$evaluate(FALSE), stats::dlnorm(y, 0, 1, FALSE))