refactor recruitment deviations to log space

inst/include/common/model.hpp

@@ -83,7 +83,7 @@ class Model {  // may need singleton
     vector<vector<Type> > naa(n_pops);
     vector<vector<Type> > ssb(n_pops);
     vector<vector<Type> > biomass(n_pops);
-    vector<vector<Type> > rec_dev(n_pops);
+    vector<vector<Type> > log_recruit_dev(n_pops);
     vector<vector<Type> > recruitment(n_pops);
     vector<vector<Type> > M(n_pops);
 #endif
@@ -148,8 +148,8 @@ class Model {  // may need singleton
 #ifdef TMB_MODEL
       naa(pop_idx) = vector<Type>((*it).second->numbers_at_age);
       ssb(pop_idx) = vector<Type>((*it).second->spawning_biomass);
-      rec_dev(pop_idx) =
-          vector<Type>((*it).second->recruitment->recruit_deviations);
+      log_recruit_dev(pop_idx) =
+          vector<Type>((*it).second->recruitment->log_recruit_devs);
       recruitment(pop_idx) = vector<Type>((*it).second->expected_recruitment);
       biomass(pop_idx) = vector<Type>((*it).second->biomass);
       M(pop_idx) = vector<Type>((*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,15 +198,15 @@ class Model {  // may need singleton
     vector<Type> NAA = ADREPORTvector(naa);
     vector<Type> Biomass = ADREPORTvector(biomass);
     vector<Type> SSB = ADREPORTvector(ssb);
-    vector<Type> RecDev = ADREPORTvector(rec_dev);
+    vector<Type> LogRecDev = ADREPORTvector(log_recruit_dev);
     vector<Type> FMort = ADREPORTvector(F_mort);
     vector<Type> ExpectedIndex = ADREPORTvector(exp_index);
     vector<Type> CNAA = ADREPORTvector(cnaa);
 
     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);

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)

inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp

@@ -28,7 +28,7 @@ class RecruitmentInterfaceBase : public FIMSRcppInterfaceBase {
   static std::map<uint32_t, RecruitmentInterfaceBase*> live_objects;
   /**< map associating the ids of RecruitmentInterfaceBase to the objects */
 
-  // static std::vector<double> recruit_deviations; /**< vector of recruitment
+  // static std::vector<double> log_recruit_devs; /**< vector of log recruitment
   // deviations*/
   // static bool constrain_deviations; /**< whether or not the rec devs are constrained*/
 
@@ -71,9 +71,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() {}
@@ -101,13 +101,13 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase {
     fims_popdy::SRBevertonHolt<double> 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();
   }
 
@@ -148,15 +148,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];
       }
     }
 

inst/include/population_dynamics/population/population.hpp

@@ -377,26 +377,30 @@ struct Population : public fims_model_object::FIMSObject<Type> {
    *
    * @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;
     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[i_dev - 1]
                    << std::endl;
     POPULATION_LOG << "rec eval = "
                    << this->recruitment->evaluate(
                           this->spawning_biomass[year - 1], phi0)
                    << std::endl;
     this->numbers_at_age[i_age_year] =
         this->recruitment->evaluate(this->spawning_biomass[year - 1], phi0) *
-        this->recruitment->recruit_deviations[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.*/
+        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 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;
   }
 
@@ -644,7 +648,7 @@ struct Population : public fims_model_object::FIMSObject<Type> {
             // 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);
 

inst/include/population_dynamics/recruitment/functors/recruitment_base.hpp

@@ -32,16 +32,16 @@ struct RecruitmentBase : public fims_model_object::FIMSObject<Type> {
   static uint32_t id_g; /**< reference id for recruitment object*/
 
   typename fims::ModelTraits<Type>::ParameterVector
-      recruit_deviations;            /**< A vector of recruitment deviations */
-  bool constrain_deviations = false; /**< A flag to indicate if recruitment
+      log_recruit_devs;            /*!< A vector of log recruitment deviations */
+  bool constrain_deviations = false; /*!< A flag to indicate if recruitment
                                  deviations are summing to zero or not */
 
   Type log_sigma_recruit; /**< Log standard deviation of log recruitment
                        deviations */
   Type log_rzero;         /**< Log of unexploited recruitment.*/
 
-  bool estimate_recruit_deviations =
-      true; /**< A flag to indicate if recruitment deviations are estimated or
+  bool estimate_log_recruit_devs =
+      true; /*!< A flag to indicate if recruitment deviations are estimated or
      not */
 
   /** @brief Constructor.
@@ -72,14 +72,14 @@ struct RecruitmentBase : public fims_model_object::FIMSObject<Type> {
   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<Type> 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> {
 
     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;
     }
   }
 };

tests/gtest/test_population_Recruitment.cpp

@@ -39,12 +39,15 @@ namespace
         // vector for storing expected recruitment 
         std::vector<double> expect_recruitment(population.nyears * population.nages, 0.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.*/
         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.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
-        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.

tests/gtest/test_population_dynamics_recruitment_base.cpp

@@ -7,16 +7,16 @@ namespace
   TEST(RecruitmentDeviations, ConstraintWorks)
   {
     fims_popdy::SRBevertonHolt<double> 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<double> 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<double> 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);
     }
   }

tests/gtest/test_population_test_fixture.hpp

@@ -142,22 +142,26 @@ class PopulationEvaluateTestFixture : public testing::Test {
     }
 
     population.growth = growth;
-
     population.Prepare();
 
+
     auto maturity = std::make_shared<fims_popdy::LogisticMaturity<double>>();
     maturity->inflection_point = 6;
     maturity->slope = 0.15;
     population.maturity = maturity;
 
-    auto recruitment = std::make_shared<fims_popdy::SRBevertonHolt<double>>();
-    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;
-    }
-    population.recruitment = recruitment;
+
+      auto recruitment = std::make_shared<fims_popdy::SRBevertonHolt<double>>();
+      recruitment->logit_steep = fims_math::logit(0.2, 1.0, 0.75);
+      recruitment->log_rzero = fims_math::log(1000000.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;
 
     int year = 4;
     int age = 6;

tests/integration/integration_class.hpp

@@ -548,18 +548,21 @@ 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);
+            /*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) {
                     JsonArray rdev = (*it).second.GetArray();
                     if (print_statements) {
                         std::cout << "recruitment deviations: ";
                     }
-                    for (size_t i = 0; i < rdev.size(); i++) {
-                        rec->recruit_deviations[i] = std::exp(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->recruit_deviations[i] << " ";
+                            std::cout << rec->log_recruit_devs[i] << " ";
                         }
                     }
                     if (print_statements) {

tests/integration/integration_test_population_tmb_nointerface.cpp

@@ -137,8 +137,11 @@ Type objective_function<Type>::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);
+  /*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;
 
   //Set maturity