From c7cc5004922eba77dfbe547abbfb22e8c9844ce7 Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA Date: Tue, 25 Jul 2023 14:28:22 -0700 Subject: [PATCH 01/11] bug: REPORT_F in model.hpp --- inst/include/common/model.hpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/inst/include/common/model.hpp b/inst/include/common/model.hpp index 72b7a23f6..e6de76672 100644 --- a/inst/include/common/model.hpp +++ b/inst/include/common/model.hpp @@ -107,6 +107,14 @@ class Model { // may need singleton jnll = rec_nll + age_comp_nll + index_nll; + + #ifdef TMB_MODEL + int i = 1; + REPORT_F(i, of); + REPORT_F(rec_nll, of); + #endif + + return jnll; } }; From 1aebd9845c3494d52be79575b050bda16d253884 Mon Sep 17 00:00:00 2001 From: timjmiller Date: Wed, 26 Jul 2023 16:02:53 -0400 Subject: [PATCH 02/11] change T to Type in model.hpp --- inst/include/common/model.hpp | 42 +++++++++++++++++------------------ 1 file changed, 21 insertions(+), 21 deletions(-) diff --git a/inst/include/common/model.hpp b/inst/include/common/model.hpp index e6de76672..09ec782bd 100644 --- a/inst/include/common/model.hpp +++ b/inst/include/common/model.hpp @@ -26,48 +26,48 @@ namespace fims { /** * @brief Model class. FIMS objective function. */ -template +template class Model { // may need singleton public: - static std::shared_ptr > + static std::shared_ptr > fims_model; /*!< Create a shared fims_model as a pointer to Model*/ - std::shared_ptr > + std::shared_ptr > fims_information; /*!< Create a shared fims_information as a pointer to Information*/ #ifdef TMB_MODEL - ::objective_function *of; + ::objective_function *of; #endif // constructor virtual ~Model() {} /** - * Returns a single Information object for type T. + * Returns a single Information object for type Type. * - * @return singleton for type T + * @return singleton for type Type */ - static std::shared_ptr > GetInstance() { - if (Model::fims_model == nullptr) { - Model::fims_model = std::make_shared >(); - Model::fims_model->fims_information = - fims::Information::GetInstance(); + static std::shared_ptr > GetInstance() { + if (Model::fims_model == nullptr) { + Model::fims_model = std::make_shared >(); + Model::fims_model->fims_information = + fims::Information::GetInstance(); } - return Model::fims_model; + return Model::fims_model; } /** * @brief Evaluate. Calculates the joint negative log-likelihood function. */ - const T Evaluate() { + const Type Evaluate() { // jnll = negative-log-likelihood (the objective function) - T jnll = 0.0; - T rec_nll = 0.0; // recrutiment nll - T age_comp_nll = 0.0; // age composition nll - T index_nll = 0.0; // survey and fishery cacth nll + Type jnll = 0.0; + Type rec_nll = 0.0; // recrutiment nll + Type age_comp_nll = 0.0; // age composition nll + Type index_nll = 0.0; // survey and fishery cacth nll // Loop over populations, evaluate, and sum up the recruitment likelihood // component - typename fims::Information::population_iterator it; + typename fims::Information::population_iterator it; for (it = this->fims_information->populations.begin(); it != this->fims_information->populations.end(); ++it) { //(*it).second points to the Population module @@ -86,7 +86,7 @@ class Model { // may need singleton FIMS_LOG << "rec nll: " << rec_nll << std::endl; } - typename fims::Information::fleet_iterator jt; + typename fims::Information::fleet_iterator jt; for (jt = this->fims_information->fleets.begin(); jt != this->fims_information->fleets.end(); ++jt) { #ifdef TMB_MODEL @@ -120,8 +120,8 @@ class Model { // may need singleton }; // Create singleton instance of Model class -template -std::shared_ptr > Model::fims_model = +template +std::shared_ptr > Model::fims_model = nullptr; // singleton instance } // namespace fims From 1259d9eeb616d7ebc0bb39a1bc8c58174399cf67 Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA Date: Thu, 7 Sep 2023 12:26:42 -0700 Subject: [PATCH 03/11] report out nll components --- inst/include/common/model.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/inst/include/common/model.hpp b/inst/include/common/model.hpp index 09ec782bd..8893eceea 100644 --- a/inst/include/common/model.hpp +++ b/inst/include/common/model.hpp @@ -109,11 +109,11 @@ class Model { // may need singleton #ifdef TMB_MODEL - int i = 1; - REPORT_F(i, of); REPORT_F(rec_nll, of); + REPORT_F(age_comp_nll, of); + REPORT_F(index_nll, of); #endif - + return jnll; } From 993e1fef22811511a1cddd1abe3bbcadb04a9317 Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA Date: Thu, 7 Sep 2023 13:19:08 -0700 Subject: [PATCH 04/11] feat: nested vector prototype for reporting out --- inst/include/common/model.hpp | 18 +++++++++++++++--- inst/include/interface/interface.hpp | 2 ++ 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/inst/include/common/model.hpp b/inst/include/common/model.hpp index 8893eceea..4e2f31ca6 100644 --- a/inst/include/common/model.hpp +++ b/inst/include/common/model.hpp @@ -65,6 +65,14 @@ class Model { // may need singleton Type rec_nll = 0.0; // recrutiment nll Type age_comp_nll = 0.0; // age composition nll Type index_nll = 0.0; // survey and fishery cacth nll + + //Create vector lists to store output for reporting + #ifdef TMB_MODEL + //VectorOfVectors creates a nested vector structure where + //each vector can be a different dimension. + typename ModelTraits::VectorOfVectors exp_index; + #endif + // Loop over populations, evaluate, and sum up the recruitment likelihood // component typename fims::Information::population_iterator it; @@ -93,16 +101,20 @@ class Model { // may need singleton (*jt).second->of = this->of; #endif age_comp_nll += (*jt).second->evaluate_age_comp_nll(); - FIMS_LOG << "survey and fleet age comp nll sum: " << age_comp_nll - << std::endl; index_nll += (*jt).second->evaluate_index_nll(); - FIMS_LOG << "survey and fleet index nll sum: " << index_nll << std::endl; + /* if ((*jt).second->is_survey == false) { #ifdef TMB_MODEL (*jt).second->of = this->of; #endif (*jt).second->ReportFleet(); } + */ + //This will not work with the current map/iterator structure but will work if + //fleets in information is defined as a vector of pointers instead and jt is an int + #ifdef TMB_MODEL + // exp_index(jt) = (*jt).second->expected_index; + #endif } jnll = rec_nll + age_comp_nll + index_nll; diff --git a/inst/include/interface/interface.hpp b/inst/include/interface/interface.hpp index 83cf46810..ba7db9deb 100644 --- a/inst/include/interface/interface.hpp +++ b/inst/include/interface/interface.hpp @@ -48,6 +48,8 @@ struct ModelTraits { of the data that is differentiable */ typedef typename CppAD::vector ParameterVector; /**< This is a vector of the parameters that is differentiable */ + typedef typename CppAD::vector> VectorOfVectors; /**< This is a vector of vectors + that becomes a list of vectors when output to R and contains values that are differentiable */ typedef typename tmbutils::vector EigenVector; /**< This is a vector as defined in TMB's namespace Eigen */ typedef typename tmbutils::matrix From 2095feb13861372495fd3d78a8e30cb5fec8d2cf Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA Date: Tue, 26 Sep 2023 13:15:07 -0700 Subject: [PATCH 05/11] move REPORT and ADREPORT to model.hpp --- inst/include/common/model.hpp | 79 +++++++++++++++---- inst/include/interface/interface.hpp | 2 - .../population_dynamics/fleet/fleet.hpp | 22 +----- .../population/population.hpp | 50 ------------ 4 files changed, 66 insertions(+), 87 deletions(-) diff --git a/inst/include/common/model.hpp b/inst/include/common/model.hpp index 4e2f31ca6..150d82eec 100644 --- a/inst/include/common/model.hpp +++ b/inst/include/common/model.hpp @@ -66,11 +66,25 @@ class Model { // may need singleton Type age_comp_nll = 0.0; // age composition nll Type index_nll = 0.0; // survey and fishery cacth nll + int n_fleets = fims_information->fleets.size(); + int n_pops = fims_information->populations.size(); + //Create vector lists to store output for reporting #ifdef TMB_MODEL - //VectorOfVectors creates a nested vector structure where + //vector< vector > creates a nested vector structure where //each vector can be a different dimension. - typename ModelTraits::VectorOfVectors exp_index; + //fleets + vector< vector > exp_index(n_fleets); + vector< vector > exp_catch(n_fleets); + vector< vector > cnaa(n_fleets); + vector< vector > cwaa(n_fleets); + vector< vector > F_mort(n_fleets); + //populations + vector< vector > naa(n_pops); + vector< vector > ssb(n_pops); + vector< vector > biomass(n_pops); + vector< vector > rec_dev(n_pops); + vector< vector > M(n_pops); #endif // Loop over populations, evaluate, and sum up the recruitment likelihood @@ -78,7 +92,7 @@ class Model { // may need singleton typename fims::Information::population_iterator it; for (it = this->fims_information->populations.begin(); it != this->fims_information->populations.end(); ++it) { - //(*it).second points to the Population module + //(*it).second points to each individual Population module FIMS_LOG << "inside pop loop" << std::endl; // Prepare recruitment (*it).second->recruitment->Prepare(); @@ -94,36 +108,73 @@ class Model { // may need singleton FIMS_LOG << "rec nll: " << rec_nll << std::endl; } + //Loop over fleets/surveys, and sum up age comp and index nlls typename fims::Information::fleet_iterator jt; for (jt = this->fims_information->fleets.begin(); jt != this->fims_information->fleets.end(); ++jt) { + + //(*jt).second points to each individual Fleet module #ifdef TMB_MODEL (*jt).second->of = this->of; #endif age_comp_nll += (*jt).second->evaluate_age_comp_nll(); index_nll += (*jt).second->evaluate_index_nll(); - /* - if ((*jt).second->is_survey == false) { -#ifdef TMB_MODEL - (*jt).second->of = this->of; -#endif - (*jt).second->ReportFleet(); - } - */ - //This will not work with the current map/iterator structure but will work if - //fleets in information is defined as a vector of pointers instead and jt is an int + } + + //Loop over populations and fleets/surveys and fill in reporting + + //initiate population index for structuring report out objects + int pop_idx = 0; + for (it = this->fims_information->populations.begin(); + it != this->fims_information->populations.end(); ++it) { + #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); + biomass(pop_idx) = vector((*it).second->biomass); + M(pop_idx) = vector((*it).second->M); + #endif + pop_idx += 1; + + } + + //initiate fleet index for structuring report out objects + int fleet_idx = 0; + for (jt = this->fims_information->fleets.begin(); + jt != this->fims_information->fleets.end(); ++jt) { #ifdef TMB_MODEL - // exp_index(jt) = (*jt).second->expected_index; + exp_index(fleet_idx) = vector((*jt).second->expected_index); + exp_catch(fleet_idx) = vector((*jt).second->expected_catch); + F_mort(fleet_idx) = vector((*jt).second->Fmort); + cnaa(fleet_idx) = vector((*jt).second->catch_numbers_at_age); + cwaa(fleet_idx) = vector((*jt).second->catch_weight_at_age); #endif + fleet_idx += 1; } jnll = rec_nll + age_comp_nll + index_nll; + //Reporting #ifdef TMB_MODEL REPORT_F(rec_nll, of); REPORT_F(age_comp_nll, of); REPORT_F(index_nll, of); + REPORT_F(jnll, of); + REPORT_F(naa, of); + ADREPORT_F(naa, of); + REPORT_F(ssb, of); + ADREPORT_F(ssb, of); + REPORT_F(rec_dev, of); + ADREPORT_F(rec_dev, of); + REPORT_F(biomass, of); + REPORT_F(M, of); + ADREPORT_F(M, of); + REPORT_F(exp_index, of); + REPORT_F(exp_catch, of); + REPORT_F(F_mort, of); + REPORT_F(cnaa, of); + REPORT_F(cwaa, of); #endif diff --git a/inst/include/interface/interface.hpp b/inst/include/interface/interface.hpp index ba7db9deb..83cf46810 100644 --- a/inst/include/interface/interface.hpp +++ b/inst/include/interface/interface.hpp @@ -48,8 +48,6 @@ struct ModelTraits { of the data that is differentiable */ typedef typename CppAD::vector ParameterVector; /**< This is a vector of the parameters that is differentiable */ - typedef typename CppAD::vector> VectorOfVectors; /**< This is a vector of vectors - that becomes a list of vectors when output to R and contains values that are differentiable */ typedef typename tmbutils::vector EigenVector; /**< This is a vector as defined in TMB's namespace Eigen */ typedef typename tmbutils::matrix diff --git a/inst/include/population_dynamics/fleet/fleet.hpp b/inst/include/population_dynamics/fleet/fleet.hpp index 16ae2d649..743d46b31 100644 --- a/inst/include/population_dynamics/fleet/fleet.hpp +++ b/inst/include/population_dynamics/fleet/fleet.hpp @@ -141,22 +141,6 @@ struct Fleet : public FIMSObject { this->Fmort[year] = fims::exp(this->log_Fmort[year]); } } - /** - * @brief Method to report out the fleet-specific derived quantities. - * - */ - void ReportFleet() { -#ifdef TMB_MODEL - // EigenVector declares a vector type from the Eigen library, which is the - // expected type for TMB's dmultinom - typename ModelTraits::EigenVector exp_index = expected_index; - REPORT_F(exp_index, of); - typename ModelTraits::EigenVector exp_catch = expected_catch; - REPORT_F(exp_catch, of); - typename ModelTraits::EigenVector F_mort = Fmort; - REPORT_F(F_mort, of); -#endif - } virtual const Type evaluate_age_comp_nll() { Type nll = 0.0; /*!< The negative log likelihood value */ @@ -200,10 +184,7 @@ struct Fleet : public FIMSObject { nll -= dmultinom.evaluate(true); } } - // typename ModelTraits::EigenVector acomp_nll; - // acomp_nll[0] = nll; - // REPORT_F(acomp_nll, of); - fims_log::get("fleet.log") << " agecomp nll: " << nll << std::endl; + #endif return nll; } @@ -226,7 +207,6 @@ struct Fleet : public FIMSObject { fims_log::get("fleet.log") << " log obs error is: " << this->log_obs_error << std::endl; fims_log::get("fleet.log") << " sd is: " << dnorm.sd << std::endl; - fims_log::get("fleet.log") << " index nll: " << nll << std::endl; #endif return nll; } diff --git a/inst/include/population_dynamics/population/population.hpp b/inst/include/population_dynamics/population/population.hpp index ca7496797..fa2021e36 100644 --- a/inst/include/population_dynamics/population/population.hpp +++ b/inst/include/population_dynamics/population/population.hpp @@ -682,56 +682,6 @@ struct Population : public FIMSObject { FIMS_LOG << "\n"; } } - // make an intermediate value in order to set multiple members (of - -#ifdef TMB_MODEL - /*Report output*/ - // REPORT_F(int(this->nages), of); //REPORT error: call of overloaded is - // ambiguous REPORT_F(int(this->nyears), of); REPORT_F(int(this->nfleets), - // of); REPORT_F(this->numbers_at_age, of); - typename ModelTraits::EigenVector naa = this->numbers_at_age; - typename ModelTraits::EigenMatrix cnaa(this->nyears * this->nages, - this->nfleets); - typename ModelTraits::EigenMatrix cwaa(this->nyears * this->nages, - this->nfleets); - typename ModelTraits::EigenVector ssb = this->spawning_biomass; - typename ModelTraits::EigenVector rec_dev = - this->recruitment->recruit_deviations; - typename ModelTraits::EigenMatrix expected_index(this->nyears, - this->nfleets); - typename ModelTraits::EigenVector recruitment = - this->expected_recruitment; - typename ModelTraits::EigenVector biomass = this->biomass; - - for (size_t fleet_ = 0; fleet_ < this->nfleets; fleet_++) { - expected_index.col(fleet_) = typename ModelTraits::EigenVector( - fleets[fleet_]->expected_index); - cnaa.col(fleet_) = typename ModelTraits::EigenVector( - fleets[fleet_]->catch_numbers_at_age); - cwaa.col(fleet_) = typename ModelTraits::EigenVector( - fleets[fleet_]->catch_weight_at_age); - } - - REPORT_F(rec_dev, of); - ADREPORT_F(rec_dev, of); - REPORT_F(naa, of); - ADREPORT_F(naa, of); - REPORT_F(cnaa, of); - REPORT_F(cwaa, of); - REPORT_F(ssb, of); - ADREPORT_F(ssb, of); - REPORT_F(expected_index, of); - REPORT_F(recruitment, of); - REPORT_F(biomass, of); - - // ADREPORT_F(this->recruitment->rzero, of); - // ADREPORT_F(this->recruitment->steep, of); can't access steep b/c not in - // recruitment_base ADREPORT_F(this->recruitment->log_sigma_recruit, of); - // ADREPORT_F(this->M, of); - // ADREPORT_F(this->maturity->slope, of); can't access slope b/c not in - // maturity base ADREPORT_F(this->maturity->median, of); can't access median - // b/c not in maturity base -#endif } }; template From 69f772bc27a1f3afa4f47222591dfbda6d47423a Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA Date: Thu, 28 Sep 2023 15:13:43 -0700 Subject: [PATCH 06/11] add recruitment and remove ADREPORT --- inst/include/common/model.hpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/inst/include/common/model.hpp b/inst/include/common/model.hpp index 150d82eec..915c4dc63 100644 --- a/inst/include/common/model.hpp +++ b/inst/include/common/model.hpp @@ -84,6 +84,7 @@ class Model { // may need singleton vector< vector > ssb(n_pops); vector< vector > biomass(n_pops); vector< vector > rec_dev(n_pops); + vector< vector > recruitment(n_pops); vector< vector > M(n_pops); #endif @@ -131,6 +132,7 @@ class Model { // may need singleton 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); + recruitment(pop_idx) = vector((*it).second->expected_recruitment); biomass(pop_idx) = vector((*it).second->biomass); M(pop_idx) = vector((*it).second->M); #endif @@ -162,14 +164,11 @@ class Model { // may need singleton REPORT_F(index_nll, of); REPORT_F(jnll, of); REPORT_F(naa, of); - ADREPORT_F(naa, of); REPORT_F(ssb, of); - ADREPORT_F(ssb, of); REPORT_F(rec_dev, of); - ADREPORT_F(rec_dev, of); + REPORT_F(recruitment, of); REPORT_F(biomass, of); REPORT_F(M, of); - ADREPORT_F(M, of); REPORT_F(exp_index, of); REPORT_F(exp_catch, of); REPORT_F(F_mort, of); From 800897805d6659ed2f5e449e6b2d61548da23c35 Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA Date: Fri, 29 Sep 2023 14:29:07 -0700 Subject: [PATCH 07/11] add ADREPORT --- inst/include/common/model.hpp | 22 +++++++++++++++++++++- inst/include/interface/interface.hpp | 19 +++++++++++++++++++ 2 files changed, 40 insertions(+), 1 deletion(-) diff --git a/inst/include/common/model.hpp b/inst/include/common/model.hpp index 915c4dc63..da86f61b6 100644 --- a/inst/include/common/model.hpp +++ b/inst/include/common/model.hpp @@ -72,7 +72,7 @@ class Model { // may need singleton //Create vector lists to store output for reporting #ifdef TMB_MODEL //vector< vector > creates a nested vector structure where - //each vector can be a different dimension. + //each vector can be a different dimension. Does not work with ADREPORT //fleets vector< vector > exp_index(n_fleets); vector< vector > exp_catch(n_fleets); @@ -174,6 +174,26 @@ class Model { // may need singleton REPORT_F(F_mort, of); REPORT_F(cnaa, of); REPORT_F(cwaa, of); + + /*ADREPORT using ADREPORTvector defined in + * inst/include/interface/interface.hpp: + * function collapses the nested vector into a single vector + */ + vector NAA = ADREPORTvector(naa); + vector Biomass = ADREPORTvector(biomass); + vector SSB = ADREPORTvector(ssb); + vector RecDev = ADREPORTvector(rec_dev); + vector FMort = ADREPORTvector(F_mort); + vector ExpectedIndex = ADREPORTvector(exp_index); + vector CNAA = ADREPORTvector(cnaa); + + ADREPORT_F(NAA, of); + ADREPORT_F(Biomass, of); + ADREPORT_F(SSB, of); + ADREPORT_F(RecDev, of); + ADREPORT_F(FMort, of); + ADREPORT_F(ExpectedIndex, of); + ADREPORT_F(CNAA, of); #endif diff --git a/inst/include/interface/interface.hpp b/inst/include/interface/interface.hpp index 83cf46810..36f3fbe95 100644 --- a/inst/include/interface/interface.hpp +++ b/inst/include/interface/interface.hpp @@ -31,6 +31,25 @@ } #define ADREPORT_F(name, F) F->reportvector.push(name, #name); +template +vector ADREPORTvector(vector< vector > x){ + int outer_dim = x.size(); + int dim = 0; + for(int i=0; i res(dim); + int idx = 0; + for(int i=0; i::value && F->do_simulate) #endif /* TMB_MODEL */ From d19ab8e17215bad82bacaf288bb7c5e4442f7953 Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA Date: Fri, 29 Sep 2023 14:30:05 -0700 Subject: [PATCH 08/11] fix tests after reporting refactor --- tests/testthat/test-fims-estimation.R | 174 ++++++++++++++++---------- 1 file changed, 111 insertions(+), 63 deletions(-) diff --git a/tests/testthat/test-fims-estimation.R b/tests/testthat/test-fims-estimation.R index f17b682fd..cb62422e9 100644 --- a/tests/testthat/test-fims-estimation.R +++ b/tests/testthat/test-fims-estimation.R @@ -167,7 +167,7 @@ test_that("deterministic test of fims", { sdr <- TMB::sdreport(obj) sdr_fixed <- summary(sdr, "fixed") - report <- obj$report() + report <- obj$report(obj$par) # log(R0) fims_logR0 <- sdr_fixed[1, "Estimate"] @@ -179,27 +179,27 @@ test_that("deterministic test of fims", { # Numbers at age for (i in 1:length(c(t(om_output$N.age)))) { - naa_are <- abs(report$naa[i] - c(t(om_output$N.age))[i]) / + naa_are <- abs(report$naa[[1]][i] - c(t(om_output$N.age))[i]) / c(t(om_output$N.age))[i] expect_lte(naa_are, 0.001) } # Biomass for (i in 1:length(om_output$biomass.mt)) { - biomass_are <- abs(report$biomass[i] - om_output$biomass.mt[i]) / + biomass_are <- abs(report$biomass[[1]][i] - om_output$biomass.mt[i]) / om_output$biomass.mt[i] expect_lte(biomass_are, 0.001) } # Spawning biomass for (i in 1:length(om_output$SSB)) { - sb_are <- abs(report$ssb[i] - om_output$SSB[i]) / + sb_are <- abs(report$ssb[[1]][i] - om_output$SSB[i]) / om_output$SSB[i] expect_lte(sb_are, 0.001) } # Recruitment - fims_naa <- matrix(report$naa[1:(om_input$nyr * om_input$nages)], + fims_naa <- matrix(report$naa[[1]][1:(om_input$nyr * om_input$nages)], nrow = om_input$nyr, byrow = TRUE ) @@ -211,46 +211,46 @@ test_that("deterministic test of fims", { expect_equal( fims_naa[1:om_input$nyr, 1], - report$recruitment[1:om_input$nyr] + report$recruitment[[1]][1:om_input$nyr] ) for (i in 1:length(om_output$N.age[, 1])) { - recruitment_are <- abs(report$recruitment[i] - om_output$N.age[i, 1]) / + recruitment_are <- abs(report$recruitment[[1]][i] - om_output$N.age[i, 1]) / om_output$N.age[i, 1] expect_lte(recruitment_are, 0.001) } # recruitment deviations (fixed at initial "true" values) - expect_equal(log(report$rec_dev), om_input$logR.resid) + expect_equal(log(report$rec_dev[[1]]), om_input$logR.resid) # F (fixed at initial "true" values) - expect_equal(report$F_mort, om_output$f) + expect_equal(report$F_mort[[1]], om_output$f) # Expected catch - fims_object <- report$expected_index[, 1] + fims_index <- report$exp_index # Expect small relative error for deterministic test for (i in 1:length(om_output$L.mt$fleet1)) { - fims_object_are <- abs(fims_object[i] - om_output$L.mt$fleet1[i]) / + fims_object_are <- abs(fims_index[[1]][i] - om_output$L.mt$fleet1[i]) / om_output$L.mt$fleet1[i] expect_lte(fims_object_are, 0.001) } 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_object[i] - em_input$L.obs$fleet1[i]) / em_input$L.obs$fleet1[i] + 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) # Expected catch number at age for (i in 1:length(c(t(om_output$L.age$fleet1)))) { - cnaa_are <- abs(report$cnaa[i, 1] - c(t(om_output$L.age$fleet1))[i]) / + cnaa_are <- abs(report$cnaa[[1]][i] - c(t(om_output$L.age$fleet1))[i]) / c(t(om_output$L.age$fleet1))[i] expect_lte(cnaa_are, 0.001) } # Expected catch number at age in proportion - fims_cnaa <- matrix(report$cnaa[1:(om_input$nyr * om_input$nages), 1], + 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) @@ -261,31 +261,29 @@ test_that("deterministic test of fims", { } # Expected survey index - fims_object <- report$expected_index[, 2] - - cwaa <- matrix(report$cwaa[1:(om_input$nyr * om_input$nages), 2], nrow = om_input$nyr, byrow = T) - expect_equal(fims_object, apply(cwaa, 1, sum) * om_output$survey_q$survey1) + cwaa <- matrix(report$cwaa[[2]][1:(om_input$nyr * om_input$nages)], nrow = om_input$nyr, byrow = T) + 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)) { - fims_object_are <- abs(fims_object[i] - om_output$survey_index_biomass$survey1[i]) / + fims_object_are <- abs(fims_index[[2]][i] - om_output$survey_index_biomass$survey1[i]) / om_output$survey_index_biomass$survey1[i] expect_lte(fims_object_are, 0.001) } 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_object[i] - em_input$surveyB.obs$survey1[i]) / em_input$surveyB.obs$survey1[i] + 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[1:(om_input$nyr * om_input$nages), 2], + 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)))) { - cnaa_are <- abs(report$cnaa[i, 2] - c(t(om_output$survey_age_comp$survey1))[i]) / + cnaa_are <- abs(report$cnaa[[2]][i] - c(t(om_output$survey_age_comp$survey1))[i]) / c(t(om_output$survey_age_comp$survey1))[i] expect_lte(cnaa_are, 0.001) } @@ -322,7 +320,7 @@ test_that("nll test of fims", { fims_logR0 <- sdr_fixed[1, "Estimate"] expect_lte(abs(fims_logR0 - log(om_input$R0)) / log(om_input$R0), 0.0001) - report <- obj$report() + report <- obj$report(obj$par) obj <- TMB::MakeADFun(data = list(), parameters, DLL = "FIMS", map = map) jnll <- obj$fn() @@ -366,6 +364,9 @@ test_that("nll test of fims", { 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() @@ -389,103 +390,134 @@ test_that("estimation test of fims", { control = list(maxit = 1000000, reltol = 1e-15) )) - report <- obj$report() + 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"), ] + 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 > 2 * sdr_naa[1:length(c(t(om_output$N.age))), 2]), + 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 (needs to be updated when std.error is available) - # sdr_biomass <- sdr_report[which(rownames(sdr_report) == "biomass"),] + # 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(report$biomass[i] - om_output$biomass.mt[i]) / om_output$biomass.mt[i] - expect_lte(biomass_are[i], 0.15) + 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 > 2 * sdr_biomass[1:length(om_output$biomass.mt),2]), - # 0.05*length(om_output$biomass.mt) - # ) + 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"), ] + 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(report$ssb[i] - om_output$SSB[i]) / om_output$SSB[i] + 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 > 2 * sdr_sb[1:length(om_output$SSB), 2]), + 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:(om_input$nyr * om_input$nages)], + 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_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, length(om_output$N.age[, 1])) for (i in 1:length(om_output$N.age[, 1])) { - 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) + 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 > 2 * sdr_naa1[1:length(om_output$SSB)]), + sum(fims_naa1_are > qnorm(.975) * sdr_naa1[1:length(om_output$SSB)]), 0.05 * length(om_output$SSB) ) expect_equal( - fims_naa[(1:om_input$nyr), 1], - report$recruitment[1:om_input$nyr] + fims_naa[, 1], + report$recruitment[[1]][1:om_input$nyr] ) # recruitment deviations - sdr_rdev <- sdr_report[which(rownames(sdr_report) == "rec_dev"), ] + sdr_rdev <- sdr_report[which(rownames(sdr_report) == "RecDev"), ] 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[i] - exp(om_input$logR.resid[i])) / - exp(om_input$logR.resid[i]) + rdev_are[i] <- abs(report$rec_dev[[1]][i] - exp(om_input$logR.resid[i]))# / + # exp(om_input$logR.resid[i]) # expect_lte(rdev_are[i], 1) # 1 } expect_lte( - sum(rdev_are > 2 * sdr_rdev[1:length(om_input$logR.resid), 2]), + sum(rdev_are > qnorm(.975) * sdr_rdev[1:length(om_input$logR.resid), 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(report$F_mort[i] - om_output$f[i]) / om_output$f[i] - expect_lte(f_are[i], 0.4) + 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 catch (needs to be updated when std.error is available) - fims_object <- report$expected_index[, 1] + # 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)) { - expect_lte(abs(fims_object[i] - om_output$L.mt$fleet1[i]) / om_output$L.mt$fleet1[i], 0.01) + 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)) { - expect_lte(abs(fims_object[i] - em_input$L.obs$fleet1[i]) / em_input$L.obs$fleet1[i], 0.025) + 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 catch number at age (needs to be updated when std.error is available) + # 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)))) { - expect_lte(abs(report$cnaa[i, 1] - c(t(om_output$L.age$fleet1))[i]) / - c(t(om_output$L.age$fleet1))[i], 0.25) + 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], @@ -503,14 +535,30 @@ test_that("estimation test of fims", { # } - # Expected survey index (needs to be updated when std.error is available) - fims_object <- report$expected_index[, 2] + # 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)) { - expect_lte(abs(fims_object[i] - om_output$survey_index_biomass$survey1[i]) / - om_output$survey_index_biomass$survey1[i], 0.05) + 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_object[i] - em_input$surveyB.obs$survey1[i]) / + expect_lte(abs(fims_survey[i,1] - em_input$surveyB.obs$survey1[i]) / em_input$surveyB.obs$survey1[i], 0.25) } From c2f1f09ccb2ef7fc2622079df4d21b70b38a2b0a Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA Date: Fri, 29 Sep 2023 15:15:19 -0700 Subject: [PATCH 09/11] fix vignette after reporting refactor --- vignettes/fims-demo.Rmd | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/vignettes/fims-demo.Rmd b/vignettes/fims-demo.Rmd index 6b3f18634..5dde51e8e 100644 --- a/vignettes/fims-demo.Rmd +++ b/vignettes/fims-demo.Rmd @@ -358,7 +358,7 @@ print(sdr_fixed) library(ggplot2) index_results <- data.frame( observed = survey_fleet_index$index_data, - expected = report$expected_index[, 2] + expected = report$expected_index[[2]] ) print(index_results) @@ -371,7 +371,7 @@ ggplot(index_results, aes(x = 1:nyears, y = observed)) + catch_results <- data.frame( observed = fishing_fleet_index$index_data, - expected = report$expected_index[, 1] + expected = report$expected_index[[1]] ) print(catch_results) From df6403a5ab81d601de6696695c6b9241634835bd Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA Date: Fri, 29 Sep 2023 15:46:47 -0700 Subject: [PATCH 10/11] fix report name --- vignettes/fims-demo.Rmd | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/vignettes/fims-demo.Rmd b/vignettes/fims-demo.Rmd index 5dde51e8e..4657a8da8 100644 --- a/vignettes/fims-demo.Rmd +++ b/vignettes/fims-demo.Rmd @@ -358,7 +358,7 @@ print(sdr_fixed) library(ggplot2) index_results <- data.frame( observed = survey_fleet_index$index_data, - expected = report$expected_index[[2]] + expected = report$exp_index[[2]] ) print(index_results) @@ -371,7 +371,7 @@ ggplot(index_results, aes(x = 1:nyears, y = observed)) + catch_results <- data.frame( observed = fishing_fleet_index$index_data, - expected = report$expected_index[[1]] + expected = report$exp_index[[1]] ) print(catch_results) From e908d9af86f1f0ac74dba3d5f746a53c6f37f84a Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA Date: Wed, 18 Oct 2023 13:54:35 -0700 Subject: [PATCH 11/11] documentation for obj$report parameter list --- tests/testthat/test-fims-estimation.R | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/tests/testthat/test-fims-estimation.R b/tests/testthat/test-fims-estimation.R index cb62422e9..8510dc3e1 100644 --- a/tests/testthat/test-fims-estimation.R +++ b/tests/testthat/test-fims-estimation.R @@ -167,6 +167,8 @@ test_that("deterministic test of fims", { 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) # log(R0) @@ -320,6 +322,8 @@ test_that("nll test of fims", { fims_logR0 <- sdr_fixed[1, "Estimate"] expect_lte(abs(fims_logR0 - log(om_input$R0)) / log(om_input$R0), 0.0001) + # 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() @@ -390,6 +394,8 @@ test_that("estimation test of fims", { 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")