rcpp_population.hpp

Go to the documentation of this file.

 
#ifndef FIMS_INTERFACE_RCPP_RCPP_OBJECTS_RCPP_POPULATION_HPP
#define FIMS_INTERFACE_RCPP_RCPP_OBJECTS_RCPP_POPULATION_HPP
 
#include "rcpp_interface_base.hpp"
#include "../../../population_dynamics/population/population.hpp"
 
class PopulationInterfaceBase : public FIMSRcppInterfaceBase {
 public:
  static uint32_t id_g;
  uint32_t id;
  static std::map<uint32_t, std::shared_ptr<PopulationInterfaceBase>>
      live_objects;
 
  SharedBoolean initialize_catch_at_age;
  SharedBoolean initialize_surplus_production;
  PopulationInterfaceBase() {
    this->id = PopulationInterfaceBase::id_g++;
    /* Create instance of map: key is id and value is pointer to
    PopulationInterfaceBase */
    // PopulationInterfaceBase::live_objects[this->id] = this;
  }
 
  PopulationInterfaceBase(const PopulationInterfaceBase &other)
      : id(other.id) {}
 
  virtual ~PopulationInterfaceBase() {}
 
  virtual uint32_t get_id() = 0;
};
 
class PopulationInterface : public PopulationInterfaceBase {
 public:
  SharedInt n_ages = 0;
  SharedInt n_fleets;
  std::shared_ptr<std::set<uint32_t>> fleet_ids;
  typedef typename std::set<uint32_t>::iterator fleet_ids_iterator;
  SharedInt n_years;
  SharedInt n_lengths;
  SharedInt maturity_id;
  SharedInt growth_id;
  SharedInt recruitment_id;
  SharedInt recruitment_err_id;
  VariableVector log_M;
  VariableVector spawning_biomass_ratio;
  VariableVector log_f_multiplier;
  VariableVector log_init_naa;
  VariableVector proportion_female;
  RealVector ages;
  SharedString name = fims::to_string("NA");
 
  // Population based derived quantities
  VariableVector total_landings_weight;
 
  VariableVector total_landings_numbers;
 
  VariableVector mortality_F;
 
  VariableVector mortality_M;
 
  VariableVector mortality_Z;
 
  VariableVector numbers_at_age;
 
  VariableVector unfished_numbers_at_age;
 
  VariableVector biomass;
  VariableVector spawning_biomass;
  VariableVector unfished_biomass;
  VariableVector unfished_spawning_biomass;
  VariableVector proportion_mature_at_age;
  VariableVector expected_recruitment;
 
  VariableVector sum_selectivity;
 
  PopulationInterface() : PopulationInterfaceBase() {
    this->proportion_female[0].initial_value_m = static_cast<double>(0.5);
    this->proportion_female[0].estimation_type_m.set("constant");
    this->fleet_ids = std::make_shared<std::set<uint32_t>>();
    std::shared_ptr<PopulationInterface> population =
        std::make_shared<PopulationInterface>(*this);
    FIMSRcppInterfaceBase::fims_interface_objects.push_back(population);
    PopulationInterfaceBase::live_objects[this->id] = population;
  }
 
  PopulationInterface(const PopulationInterface &other)
      : PopulationInterfaceBase(other),
        n_ages(other.n_ages),
        n_fleets(other.n_fleets),
        fleet_ids(other.fleet_ids),
        n_years(other.n_years),
        n_lengths(other.n_lengths),
        maturity_id(other.maturity_id),
        growth_id(other.growth_id),
        recruitment_id(other.recruitment_id),
        recruitment_err_id(other.recruitment_id),
        log_M(other.log_M),
        spawning_biomass_ratio(other.spawning_biomass_ratio),
        log_f_multiplier(other.log_f_multiplier),
        log_init_naa(other.log_init_naa),
        proportion_female(other.proportion_female),
        ages(other.ages),
        name(other.name),
        total_landings_weight(other.total_landings_weight),
        total_landings_numbers(other.total_landings_numbers),
        mortality_F(other.mortality_F),
        mortality_M(other.mortality_M),
        mortality_Z(other.mortality_Z),
        numbers_at_age(other.numbers_at_age),
        unfished_numbers_at_age(other.unfished_numbers_at_age),
        biomass(other.biomass),
        spawning_biomass(other.spawning_biomass),
        unfished_biomass(other.unfished_biomass),
        unfished_spawning_biomass(other.unfished_spawning_biomass),
        proportion_mature_at_age(other.proportion_mature_at_age),
        expected_recruitment(other.expected_recruitment),
        sum_selectivity(other.sum_selectivity) {}
 
  virtual ~PopulationInterface() {}
 
  virtual uint32_t get_id() { return this->id; }
 
  void SetName(const std::string &name) { this->name.set(name); }
 
  std::string GetName() const { return this->name.get(); }
 
  void SetMaturityID(uint32_t maturity_id) {
    this->maturity_id.set(maturity_id);
  }
 
  void SetGrowthID(uint32_t growth_id) { this->growth_id.set(growth_id); }
 
  void SetRecruitmentID(uint32_t recruitment_id) {
    this->recruitment_id.set(recruitment_id);
  }
 
  void AddFleet(uint32_t fleet_id) { this->fleet_ids->insert(fleet_id); }
 
  virtual void finalize() {
    if (this->finalized) {
      // log warning that finalize has been called more than once.
      FIMS_WARNING_LOG("Population " + fims::to_string(this->id) +
                       " has been finalized already.");
    }
 
    this->finalized = true;  // indicate this has been called already
 
    std::shared_ptr<fims_info::Information<double>> info =
        fims_info::Information<double>::GetInstance();
 
    fims_info::Information<double>::population_iterator it;
 
    it = info->populations.find(this->id);
 
    std::shared_ptr<fims_popdy::Population<double>> pop =
        info->populations[this->id];
    it = info->populations.find(this->id);
    if (it == info->populations.end()) {
      FIMS_WARNING_LOG("Population " + fims::to_string(this->id) +
                       " not found in Information.");
      return;
    } else {
      for (size_t i = 0; i < this->log_M.size(); i++) {
        if (this->log_M[i].estimation_type_m.get() == "constant") {
          this->log_M[i].final_value_m = this->log_M[i].initial_value_m;
        } else {
          this->log_M[i].final_value_m = pop->log_M[i];
        }
      }
 
      for (size_t i = 0; i < this->log_f_multiplier.size(); i++) {
        if (this->log_f_multiplier[i].estimation_type_m.get() == "constant") {
          this->log_f_multiplier[i].final_value_m =
              this->log_f_multiplier[i].initial_value_m;
        } else {
          this->log_f_multiplier[i].final_value_m = pop->log_f_multiplier[i];
        }
      }
 
      for (size_t i = 0; i < this->log_init_naa.size(); i++) {
        if (this->log_init_naa[i].estimation_type_m.get() == "constant") {
          this->log_init_naa[i].final_value_m =
              this->log_init_naa[i].initial_value_m;
        } else {
          this->log_init_naa[i].final_value_m = pop->log_init_naa[i];
        }
      }
 
      for (size_t i = 0; i < this->proportion_female.size(); i++) {
        if (this->proportion_female[i].estimation_type_m.get() == "constant") {
          this->proportion_female[i].final_value_m =
              this->proportion_female[i].initial_value_m;
        } else {
          this->proportion_female[i].final_value_m =
              pop->proportion_female.get_force_scalar(i);
        }
      }
    }
  }
 
#ifdef TMB_MODEL
 
  template <typename Type>
  bool add_to_fims_tmb_internal() {
    std::shared_ptr<fims_info::Information<Type>> info =
        fims_info::Information<Type>::GetInstance();
 
    std::shared_ptr<fims_popdy::Population<Type>> population =
        std::make_shared<fims_popdy::Population<Type>>();
 
    std::stringstream ss;
 
    // set relative info
    population->id = this->id;
    population->n_years = this->n_years.get();
    population->n_fleets = this->n_fleets.get();
    // only define ages if n_ages greater than 0
    if (this->n_ages.get() > 0) {
      population->n_ages = this->n_ages.get();
      if (static_cast<size_t>(this->n_ages.get()) == this->ages.size()) {
        population->ages.resize(this->n_ages.get());
      } else {
        throw std::invalid_argument(
            "The size of the ages vector for population " +
            fims::to_string(this->id) + " is not equal to n_ages.");
      }
    }
 
    fleet_ids_iterator it;
    for (it = this->fleet_ids->begin(); it != this->fleet_ids->end(); it++) {
      population->fleet_ids.insert(*it);
    }
 
    population->growth_id = this->growth_id.get();
    population->recruitment_id = this->recruitment_id.get();
    population->maturity_id = this->maturity_id.get();
    population->log_M.resize(this->log_M.size());
 
    if (this->log_f_multiplier.size() ==
        static_cast<size_t>(this->n_years.get())) {
      population->log_f_multiplier.resize(this->log_f_multiplier.size());
    } else {
      FIMS_WARNING_LOG(
          "The log_f_multiplier vector is not of size n_years. Filling with "
          "zeros.");
      this->log_f_multiplier.resize((this->n_years.get()));
      for (size_t i = 0; i < log_f_multiplier.size(); i++) {
        this->log_f_multiplier[i].initial_value_m = static_cast<double>(0.0);
        this->log_f_multiplier[i].estimation_type_m.set("constant");
      }
      population->log_f_multiplier.resize(this->log_f_multiplier.size());
    }
 
    if (this->spawning_biomass_ratio.size() ==
        static_cast<size_t>(this->n_years.get() + 1)) {
      population->spawning_biomass_ratio.resize(
          this->spawning_biomass_ratio.size());
    } else {
      FIMS_WARNING_LOG(
          "Setting spawning_biomass_ratio vector to size n_years + 1.");
      this->spawning_biomass_ratio.resize((this->n_years.get() + 1));
      population->spawning_biomass_ratio.resize(
          this->spawning_biomass_ratio.size());
    }
    info->variable_map[this->spawning_biomass_ratio.id_m] =
        &(population)->spawning_biomass_ratio;
 
    population->log_init_naa.resize(this->log_init_naa.size());
    for (size_t i = 0; i < log_M.size(); i++) {
      population->log_M[i] = this->log_M[i].initial_value_m;
      if (this->log_M[i].estimation_type_m.get() == "fixed_effects") {
        ss.str("");
        ss << "Population." << this->id << ".log_M." << this->log_M[i].id_m;
        info->RegisterParameterName(ss.str());
        info->RegisterParameter(population->log_M[i]);
      }
      if (this->log_M[i].estimation_type_m.get() == "random_effects") {
        ss.str("");
        ss << "Population." << this->id << ".log_M." << this->log_M[i].id_m;
        info->RegisterRandomEffectName(ss.str());
        info->RegisterRandomEffect(population->log_M[i]);
      }
    }
    info->variable_map[this->log_M.id_m] = &(population)->log_M;
 
    for (size_t i = 0; i < log_f_multiplier.size(); i++) {
      population->log_f_multiplier[i] =
          this->log_f_multiplier[i].initial_value_m;
      if (this->log_f_multiplier[i].estimation_type_m.get() ==
          "fixed_effects") {
        ss.str("");
        ss << "Population." << this->id << ".log_f_multiplier."
           << this->log_f_multiplier[i].id_m;
        info->RegisterParameterName(ss.str());
        info->RegisterParameter(population->log_f_multiplier[i]);
      }
      if (this->log_f_multiplier[i].estimation_type_m.get() ==
          "random_effects") {
        ss.str("");
        ss << "Population." << this->id << ".log_f_multiplier."
           << this->log_f_multiplier[i].id_m;
        info->RegisterRandomEffectName(ss.str());
        info->RegisterRandomEffect(population->log_f_multiplier[i]);
      }
    }
    info->variable_map[this->log_f_multiplier.id_m] =
        &(population)->log_f_multiplier;
 
    for (size_t i = 0; i < log_init_naa.size(); i++) {
      population->log_init_naa[i] = this->log_init_naa[i].initial_value_m;
      if (this->log_init_naa[i].estimation_type_m.get() == "fixed_effects") {
        ss.str("");
        ss << "Population." << this->id << ".log_init_naa."
           << this->log_init_naa[i].id_m;
        info->RegisterParameterName(ss.str());
        info->RegisterParameter(population->log_init_naa[i]);
      }
      if (this->log_init_naa[i].estimation_type_m.get() == "random_effects") {
        ss.str("");
        ss << "Population." << this->id << ".log_init_naa."
           << this->log_init_naa[i].id_m;
        info->RegisterRandomEffectName(ss.str());
        info->RegisterRandomEffect(population->log_init_naa[i]);
      }
    }
    info->variable_map[this->log_init_naa.id_m] = &(population)->log_init_naa;
 
    if (this->proportion_female.size() == 1 ||
        this->proportion_female.size() ==
            static_cast<size_t>(this->n_ages.get())) {
      population->proportion_female.resize(this->proportion_female.size());
    } else {
      FIMS_WARNING_LOG(
          "The proportion_female vector is not of size 1 or n_ages. Filling "
          "with 0.5.");
      this->proportion_female.resize(1);
      this->proportion_female[0].initial_value_m = static_cast<double>(0.5);
      this->proportion_female[0].estimation_type_m.set("constant");
      population->proportion_female.resize(this->proportion_female.size());
    }
 
    for (size_t i = 0; i < this->proportion_female.size(); i++) {
      if (this->proportion_female[i].initial_value_m < 0.0 ||
          this->proportion_female[i].initial_value_m > 1.0) {
        FIMS_WARNING_LOG(
            "proportion_female should be in [0, 1]; got " +
            fims::to_string(this->proportion_female[i].initial_value_m) +
            " at index " + fims::to_string(i) + ".");
      }
      population->proportion_female[i] =
          this->proportion_female[i].initial_value_m;
      if (this->proportion_female[i].estimation_type_m.get() ==
          "fixed_effects") {
        ss.str("");
        ss << "Population." << this->id << ".proportion_female."
           << this->proportion_female[i].id_m;
        info->RegisterParameterName(ss.str());
        info->RegisterParameter(population->proportion_female[i]);
      }
      if (this->proportion_female[i].estimation_type_m.get() ==
          "random_effects") {
        ss.str("");
        ss << "Population." << this->id << ".proportion_female."
           << this->proportion_female[i].id_m;
        info->RegisterRandomEffectName(ss.str());
        info->RegisterRandomEffect(population->proportion_female[i]);
      }
    }
    info->variable_map[this->proportion_female.id_m] =
        &(population)->proportion_female;
 
    for (size_t i = 0; i < ages.size(); i++) {
      population->ages[i] = this->ages[i];
    }
 
    // add to Information
    info->populations[population->id] = population;
 
    return true;
  }
 
  virtual bool add_to_fims_tmb() {
    this->add_to_fims_tmb_internal<TMB_FIMS_REAL_TYPE>();
    this->add_to_fims_tmb_internal<TMBAD_FIMS_TYPE>();
 
    return true;
  }
 
#endif
};
 
#endif