catch_at_age.hpp

Go to the documentation of this file.

#ifndef FIMS_MODELS_CATCH_AT_AGE_HPP
#define FIMS_MODELS_CATCH_AT_AGE_HPP
 
#include <algorithm>
#include <cmath>
#include <set>
#include <regex>
#include <stdexcept>
 
#include "fishery_model_base.hpp"
#include "population_dynamics/age_to_length_conversion/functors/runtime.hpp"
 
/* Dictionary block for shared parameter snippet documentations.
 * Referenced in function docs via @snippet{doc} this snippet_id.
  [param_population]
  @param population Shared pointer to the population object.
  [param_population]
  [param_i_age_year]
  @param i_age_year Dimension folded index for age and year.
  [param_i_age_year]
  [param_year]
  @param year Year index.
  [param_year]
  [param_age]
  @param age Age index.
  [param_age]
  [param_i_agem1_yearm1]
  @param i_agem1_yearm1 Dimension folded index for age-1 and year-1.
  [param_i_agem1_yearm1]
  [param_i_dev]
  @param i_dev Index to log_recruit_dev of vector length n_years-1.
  [param_i_dev]
  [param_other]
  @param other The other CatchAtAge object to copy from.
  [param_other]
 */
 
namespace fims_popdy {
 
template <typename Type>
class CatchAtAge : public FisheryModelBase<Type> {
 public:
  std::string name_m;

  typedef typename std::map<std::string, fims::Vector<Type>>::iterator
      derived_quantities_iterator;

  typedef typename std::map<uint32_t,
                            std::map<std::string, fims::Vector<Type>>>::iterator
      fleet_derived_quantities_iterator;

  typedef
      typename std::map<uint32_t,
                        std::map<std::string, fims::Vector<size_t>>>::iterator
          fleet_derived_quantities_dims_iterator;
  typedef typename std::map<uint32_t,
                            std::map<std::string, fims::Vector<Type>>>::iterator
      population_derived_quantities_iterator;

  typedef
      typename std::map<uint32_t,
                        std::map<std::string, fims::Vector<size_t>>>::iterator
          population_derived_quantities_dims_iterator;

  typedef typename std::map<uint32_t,
                            std::shared_ptr<fims_popdy::Fleet<Type>>>::iterator
      fleet_iterator;
  typedef
      typename std::map<std::string, fims::Vector<Type>>::iterator dq_iterator;
  std::map<std::string, fims::Vector<fims::Vector<Type>>> report_vectors;
  bool report_age_to_length_conversion_derived_tensor = false;
 
 public:
  std::vector<Type> ages; 
  CatchAtAge() : FisheryModelBase<Type>() {
    std::stringstream ss;
    ss << "caa_" << this->GetId() << "_";
    this->name_m = ss.str();
    this->model_type_m = "caa";
  }

  CatchAtAge(const CatchAtAge &other)
      : FisheryModelBase<Type>(other),
        name_m(other.name_m),
        ages(other.ages),
        report_age_to_length_conversion_derived_tensor(
            other.report_age_to_length_conversion_derived_tensor) {
    this->model_type_m = "caa";
  }

  virtual ~CatchAtAge() {}

  virtual void Initialize() {
    for (size_t p = 0; p < this->populations.size(); p++) {
      if (this->populations[p]->proportion_female.size() == 0) {
        this->populations[p]->proportion_female.resize(1);
        this->populations[p]->proportion_female[0] = static_cast<Type>(0.5);
      }
 
      this->populations[p]->M.resize(this->populations[p]->n_years *
                                     this->populations[p]->n_ages);
 
      this->populations[p]->f_multiplier.resize(this->populations[p]->n_years);
 
      this->populations[p]->spawning_biomass_ratio.resize(
          (this->populations[p]->n_years + 1));
    }
 
    for (fleet_iterator fit = this->fleets.begin(); fit != this->fleets.end();
         ++fit) {
      std::shared_ptr<fims_popdy::Fleet<Type>> &fleet = (*fit).second;
 
      if (fleet->log_q.size() == 0) {
        fleet->log_q.resize(1);
        fleet->log_q[0] = static_cast<Type>(0.0);
      }
      fleet->q.resize(fleet->log_q.size());
      fleet->Fmort.resize(fleet->n_years);
    }
  }

  virtual void Prepare() {
    for (size_t p = 0; p < this->populations.size(); p++) {
      std::shared_ptr<fims_popdy::Population<Type>> &population =
          this->populations[p];
 
      auto &derived_quantities =
          this->GetPopulationDerivedQuantities(population->GetId());
 
      // Reset the derived quantities for the population
      for (auto &kv : derived_quantities) {
        this->ResetVector(kv.second);
      }
 
      // Transformation Section
      for (size_t age = 0; age < population->n_ages; age++) {
        for (size_t year = 0; year < population->n_years; year++) {
          size_t i_age_year = age * population->n_years + year;
          population->M[i_age_year] =
              fims_math::exp(population->log_M[i_age_year]);
        }
      }
 
      for (size_t year = 0; year < population->n_years; year++) {
        population->f_multiplier[year] =
            fims_math::exp(population->log_f_multiplier[year]);
      }
    }
 
    for (fleet_iterator fit = this->fleets.begin(); fit != this->fleets.end();
         ++fit) {
      std::shared_ptr<fims_popdy::Fleet<Type>> &fleet = (*fit).second;
      auto &derived_quantities =
          this->GetFleetDerivedQuantities(fleet->GetId());
 
      for (auto &kv : derived_quantities) {
        this->ResetVector(kv.second);
      }
 
      // Transformation Section
      for (size_t i = 0; i < fleet->log_q.size(); i++) {
        fleet->q[i] = fims_math::exp(fleet->log_q[i]);
      }
 
      for (size_t year = 0; year < fleet->n_years; year++) {
        fleet->Fmort[year] = fims_math::exp(fleet->log_Fmort[year]);
      }
    }
  }

  void AddPopulation(uint32_t id) { this->population_ids.insert(id); }

  std::set<uint32_t> &GetPopulationIds() { return this->population_ids; }

  std::vector<std::shared_ptr<fims_popdy::Population<Type>>> &GetPopulations() {
    return this->populations;
  }

  void CalculateInitialNumbersAA(
      std::shared_ptr<fims_popdy::Population<Type>> &population,
      size_t i_age_year, size_t age) {
    std::map<std::string, fims::Vector<Type>> &dq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    dq_["numbers_at_age"][i_age_year] =
        fims_math::exp(population->log_init_naa[age]);
  }

  void CalculateNumbersAA(
      std::shared_ptr<fims_popdy::Population<Type>> &population,
      size_t i_age_year, size_t i_agem1_yearm1, size_t age) {
    // using Z from previous age/year
 
    std::map<std::string, fims::Vector<Type>> &dq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    dq_["numbers_at_age"][i_age_year] =
        dq_["numbers_at_age"][i_agem1_yearm1] *
        (fims_math::exp(-dq_["mortality_Z"][i_agem1_yearm1]));
 
    // Plus group calculation
    if (age == (population->n_ages - 1)) {
      dq_["numbers_at_age"][i_age_year] =
          dq_["numbers_at_age"][i_age_year] +
          dq_["numbers_at_age"][i_agem1_yearm1 + 1] *
              (fims_math::exp(-dq_["mortality_Z"][i_agem1_yearm1 + 1]));
    }
  }

  void CalculateUnfishedNumbersAA(
      std::shared_ptr<fims_popdy::Population<Type>> &population,
      size_t i_age_year, size_t i_agem1_yearm1, size_t age) {
    std::map<std::string, fims::Vector<Type>> &dq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    // using M from previous age/year
    dq_["unfished_numbers_at_age"][i_age_year] =
        dq_["unfished_numbers_at_age"][i_agem1_yearm1] *
        (fims_math::exp(-population->M[i_agem1_yearm1]));
 
    // Plus group calculation
    if (age == (population->n_ages - 1)) {
      dq_["unfished_numbers_at_age"][i_age_year] =
          dq_["unfished_numbers_at_age"][i_age_year] +
          dq_["unfished_numbers_at_age"][i_agem1_yearm1 + 1] *
              (fims_math::exp(-population->M[i_agem1_yearm1 + 1]));
    }
  }

  void CalculateUnfishedInitialNumbersAA(
      std::shared_ptr<fims_popdy::Population<Type>> &population,
      size_t i_age_year, size_t age) {
    std::map<std::string, fims::Vector<Type>> &dq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    // In year 0, the previous age is the previous element of the same year
    size_t i_agem1 = i_age_year - 1;
    dq_["unfished_numbers_at_age"][i_age_year] =
        dq_["unfished_numbers_at_age"][i_agem1] *
        fims_math::exp(-population->M[i_agem1]);
 
    if (age == (population->n_ages - 1)) {
      dq_["unfished_numbers_at_age"][i_age_year] =
          dq_["unfished_numbers_at_age"][i_age_year] /
          (1.0 - fims_math::exp(-population->M[i_age_year]));
    }
  }

  void CalculateMortality(
      std::shared_ptr<fims_popdy::Population<Type>> &population,
      size_t i_age_year, size_t year, size_t age) {
    std::map<std::string, fims::Vector<Type>> &dq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    for (size_t fleet_ = 0; fleet_ < population->n_fleets; fleet_++) {
      // evaluate is a member function of the selectivity class
      Type s = population->fleets[fleet_]->selectivity->evaluate(
          population->ages[age], year);
 
      dq_["mortality_F"][i_age_year] +=
          population->fleets[fleet_]->Fmort[year] *
          population->f_multiplier[year] * s;
 
      dq_["sum_selectivity"][i_age_year] += s;
    }
    dq_["mortality_M"][i_age_year] = population->M[i_age_year];
 
    dq_["mortality_Z"][i_age_year] =
        population->M[i_age_year] + dq_["mortality_F"][i_age_year];
  }

  void CalculateBiomass(
      std::shared_ptr<fims_popdy::Population<Type>> &population,
      size_t i_age_year, size_t year, size_t age) {
    std::map<std::string, fims::Vector<Type>> &dq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    dq_["biomass"][year] += dq_["numbers_at_age"][i_age_year] *
                            PopulationMeanWeightAA(population, year, age);
  }

  void CalculateUnfishedBiomass(
      std::shared_ptr<fims_popdy::Population<Type>> &population,
      size_t i_age_year, size_t year, size_t age) {
    std::map<std::string, fims::Vector<Type>> &dq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    dq_["unfished_biomass"][year] +=
        dq_["unfished_numbers_at_age"][i_age_year] *
        PopulationMeanWeightAA(population, year, age);
  }

  void CalculateSpawningBiomass(
      std::shared_ptr<fims_popdy::Population<Type>> &population,
      size_t i_age_year, size_t year, size_t age) {
    std::map<std::string, fims::Vector<Type>> &dq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    dq_["spawning_biomass"][year] +=
        population->proportion_female.get_force_scalar(age) *
        dq_["numbers_at_age"][i_age_year] *
        dq_["proportion_mature_at_age"][i_age_year] *
        PopulationMeanWeightAA(population, year, age);
  }

  void CalculateUnfishedSpawningBiomass(
      std::shared_ptr<fims_popdy::Population<Type>> &population,
      size_t i_age_year, size_t year, size_t age) {
    std::map<std::string, fims::Vector<Type>> &dq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    dq_["unfished_spawning_biomass"][year] +=
        population->proportion_female.get_force_scalar(age) *
        dq_["unfished_numbers_at_age"][i_age_year] *
        dq_["proportion_mature_at_age"][i_age_year] *
        PopulationMeanWeightAA(population, year, age);
  }

  void CalculateSpawningBiomassRatio(
      std::shared_ptr<fims_popdy::Population<Type>> &population, size_t year) {
    std::map<std::string, fims::Vector<Type>> &dq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
    population->spawning_biomass_ratio[year] =
        dq_["spawning_biomass"][year] / dq_["unfished_spawning_biomass"][0];
  }

  Type CalculateSBPR0(
      std::shared_ptr<fims_popdy::Population<Type>> &population) {
    std::map<std::string, fims::Vector<Type>> &dq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    std::vector<Type> numbers_spr(population->n_ages, 1.0);
    Type phi_0 = 0.0;
    phi_0 += numbers_spr[0] *
             population->proportion_female.get_force_scalar(0) *
             dq_["proportion_mature_at_age"][0] *
             PopulationMeanWeightAA(population, 0, 0);
    for (size_t a = 1; a < (population->n_ages - 1); a++) {
      numbers_spr[a] =
          numbers_spr[a - 1] * fims_math::exp(-population->M[a - 1]);
      phi_0 += numbers_spr[a] *
               population->proportion_female.get_force_scalar(a) *
               dq_["proportion_mature_at_age"][a] *
               PopulationMeanWeightAA(population, 0, a);
    }
 
    numbers_spr[population->n_ages - 1] =
        (numbers_spr[population->n_ages - 2] *
         fims_math::exp(-population->M[population->n_ages - 2])) /
        (1 - fims_math::exp(-population->M[population->n_ages - 1]));
    phi_0 +=
        numbers_spr[population->n_ages - 1] *
        population->proportion_female.get_force_scalar(population->n_ages - 1) *
        dq_["proportion_mature_at_age"][population->n_ages - 1] *
        PopulationMeanWeightAA(population, 0, population->n_ages - 1);
 
    return phi_0;
  }

  void CalculateRecruitment(
      std::shared_ptr<fims_popdy::Population<Type>> &population,
      size_t i_age_year, size_t year, size_t i_dev) {
    std::map<std::string, fims::Vector<Type>> &dq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    Type phi_0 = CalculateSBPR0(population);
 
    if (i_dev == population->n_years) {
      dq_["numbers_at_age"][i_age_year] =
          population->recruitment->evaluate_mean(
              dq_["spawning_biomass"][year - 1], phi_0);
      /*the final year of the time series has no data to inform recruitment
      devs, so this value is set to the mean recruitment.*/
    } else {
      // Why are we using evaluate_mean, how come a virtual function was
      // changed? AMH: there are now two virtual functions: evaluate_mean and
      // evaluate_process (see below)
      population->recruitment->log_expected_recruitment[year - 1] =
          fims_math::log(population->recruitment->evaluate_mean(
              dq_["spawning_biomass"][year - 1], phi_0));
 
      dq_["numbers_at_age"][i_age_year] = fims_math::exp(
          population->recruitment->process->evaluate_process(year - 1));
    }
 
    dq_["expected_recruitment"][year] = dq_["numbers_at_age"][i_age_year];
  }

  void CalculateMaturityAA(
      std::shared_ptr<fims_popdy::Population<Type>> &population,
      size_t i_age_year, size_t year, size_t age) {
    std::map<std::string, fims::Vector<Type>> &dq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    // Year n_years is the start-of-year state after the last data year, so it
    // reuses that year's maturity.
    dq_["proportion_mature_at_age"][i_age_year] =
        population->maturity->evaluate(
            population->ages[age], (std::min)(year, population->n_years - 1));
  }

  void CalculateCatch(std::shared_ptr<fims_popdy::Population<Type>> &population,
                      size_t year, size_t age) {
    std::map<std::string, fims::Vector<Type>> &pdq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    for (size_t fleet_ = 0; fleet_ < population->n_fleets; fleet_++) {
      std::map<std::string, fims::Vector<Type>> &fdq_ =
          this->GetFleetDerivedQuantities(population->fleets[fleet_]->GetId());
      size_t i_age_year = year * population->n_ages + age;
 
      pdq_["total_catch_weight"][year] +=
          fdq_["catch_weight_at_age"][i_age_year];
 
      fdq_["catch_weight"][year] += fdq_["catch_weight_at_age"][i_age_year];
 
      pdq_["total_catch_numbers"][year] +=
          fdq_["catch_numbers_at_age"][i_age_year];
 
      fdq_["catch_numbers"][year] += fdq_["catch_numbers_at_age"][i_age_year];
    }
  }

  void CalculateCatchWeightAA(
      std::shared_ptr<fims_popdy::Population<Type>> &population, size_t year,
      size_t age) {
    int i_age_year = year * population->n_ages + age;
 
    if (population->growth == nullptr) {
      throw std::runtime_error(
          "Population growth pointer was null while resolving landings "
          "weight-at-age.");
    }
 
    for (size_t fleet_ = 0; fleet_ < population->n_fleets; fleet_++) {
      std::shared_ptr<fims_popdy::Fleet<Type>> &fleet =
          population->fleets[fleet_];
      std::map<std::string, fims::Vector<Type>> &fdq_ =
          this->GetFleetDerivedQuantities(fleet->GetId());
 
      Type mean_weight_at_age;
      if (fleet->weight_at_age.size() > 0) {
        mean_weight_at_age = fleet->weight_at_age[i_age_year];
      } else {
        mean_weight_at_age = PopulationMeanWeightAA(population, year, age);
      }
 
      fdq_["catch_weight_at_age"][i_age_year] =
          fdq_["catch_numbers_at_age"][i_age_year] * mean_weight_at_age;
    }
  }

  void CalculateCatchNumbersAA(
      std::shared_ptr<fims_popdy::Population<Type>> &population,
      size_t i_age_year, size_t year, size_t age) {
    std::map<std::string, fims::Vector<Type>> &pdq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    for (size_t fleet_ = 0; fleet_ < population->n_fleets; fleet_++) {
      std::map<std::string, fims::Vector<Type>> &fdq_ =
          this->GetFleetDerivedQuantities(population->fleets[fleet_]->GetId());
 
      // Baranov Catch Equation
      fdq_["catch_numbers_at_age"][i_age_year] +=
          (population->fleets[fleet_]->Fmort[year] *
           population->f_multiplier[year] *
           population->fleets[fleet_]->selectivity->evaluate(
               population->ages[age], year)) /
          pdq_["mortality_Z"][i_age_year] * pdq_["numbers_at_age"][i_age_year] *
          (1 - fims_math::exp(-(pdq_["mortality_Z"][i_age_year])));
    }
  }

  void CalculateIndex(std::shared_ptr<fims_popdy::Population<Type>> &population,
                      size_t i_age_year, size_t year, size_t age) {
    for (size_t fleet_ = 0; fleet_ < population->n_fleets; fleet_++) {
      std::map<std::string, fims::Vector<Type>> &fdq_ =
          this->GetFleetDerivedQuantities(population->fleets[fleet_]->GetId());
 
      fdq_["index_weight"][year] += fdq_["index_weight_at_age"][i_age_year];
 
      fdq_["index_numbers"][year] += fdq_["index_numbers_at_age"][i_age_year];
    }
  }

  void CalculateIndexNumbersAA(
      std::shared_ptr<fims_popdy::Population<Type>> &population,
      size_t i_age_year, size_t year, size_t age) {
    std::map<std::string, fims::Vector<Type>> &pdq_ =
        this->GetPopulationDerivedQuantities(population->GetId());
 
    for (size_t fleet_ = 0; fleet_ < population->n_fleets; fleet_++) {
      std::map<std::string, fims::Vector<Type>> &fdq_ =
          this->GetFleetDerivedQuantities(population->fleets[fleet_]->GetId());
 
      fdq_["index_numbers_at_age"][i_age_year] +=
          (population->fleets[fleet_]->q.get_force_scalar(year) *
           population->fleets[fleet_]->selectivity->evaluate(
               population->ages[age], year)) *
          pdq_["numbers_at_age"][i_age_year];
    }
  }

  void CalculateIndexWeightAA(
      std::shared_ptr<fims_popdy::Population<Type>> &population, size_t year,
      size_t age) {
    int i_age_year = year * population->n_ages + age;
 
    if (population->growth == nullptr) {
      throw std::runtime_error(
          "Population growth pointer was null while resolving index "
          "weight-at-age.");
    }
 
    for (size_t fleet_ = 0; fleet_ < population->n_fleets; fleet_++) {
      std::shared_ptr<fims_popdy::Fleet<Type>> &fleet =
          population->fleets[fleet_];
      std::map<std::string, fims::Vector<Type>> &fdq_ =
          this->GetFleetDerivedQuantities(fleet->GetId());
 
      Type mean_weight_at_age;
      if (fleet->weight_at_age.size() > 0) {
        mean_weight_at_age = fleet->weight_at_age[i_age_year];
      } else {
        mean_weight_at_age = PopulationMeanWeightAA(population, year, age);
      }
 
      fdq_["index_weight_at_age"][i_age_year] =
          fdq_["index_numbers_at_age"][i_age_year] * mean_weight_at_age;
    }
  }

  void evaluate_age_comp() {
    fleet_iterator fit;
    for (fit = this->fleets.begin(); fit != this->fleets.end(); ++fit) {
      std::map<std::string, fims::Vector<Type>> &fdq_ =
          this->GetFleetDerivedQuantities((*fit).second->GetId());
 
      std::shared_ptr<fims_popdy::Fleet<Type>> &fleet = (*fit).second;
      for (size_t y = 0; y < fleet->n_years; y++) {
        Type sum = static_cast<Type>(0.0);
        Type sum_obs = static_cast<Type>(0.0);
        // robust_add is a small value to add to expected composition
        // proportions at age to stabilize likelihood calculations
        // when the expected proportions are close to zero.
        // Type robust_add = static_cast<Type>(0.0); // zeroed out before
        // testing 0.0001; sum robust is used to calculate the total sum of
        // robust additions to ensure that proportions sum to 1. Type robust_sum
        // = static_cast<Type>(1.0);
 
        for (size_t a = 0; a < fleet->n_ages; a++) {
          size_t i_age_year = y * fleet->n_ages + a;
          // Here we have a check to determine if the age comp
          // should be calculated from the retained catch or
          // the total population. These values are slightly different.
          // In the future this will have more impact as we implement
          // timing rather than everything occurring at the start of
          // the year.
          if (fleet->fleet_observed_catch_data_id_m == -999) {
            fdq_["agecomp_expected"][i_age_year] =
                fdq_["index_numbers_at_age"][i_age_year];
          } else {
            fdq_["agecomp_expected"][i_age_year] =
                fdq_["catch_numbers_at_age"][i_age_year];
          }
          sum += fdq_["agecomp_expected"][i_age_year];
          // robust_sum -= robust_add;
 
          // This sums over the observed age composition data so that
          // the expected age composition can be rescaled to match the
          // total number observed. The check for na values should not
          // be needed as individual years should not have missing data.
          // This is need to be re-explored if/when we modify FIMS to
          // allow for composition bins that do not match the population
          // bins.
          if (fleet->fleet_observed_agecomp_data_id_m != -999) {
            if (fleet->observed_agecomp_data->at(i_age_year) !=
                fleet->observed_agecomp_data->na_value) {
              sum_obs += fleet->observed_agecomp_data->at(i_age_year);
            }
          }
        }
        for (size_t a = 0; a < fleet->n_ages; a++) {
          size_t i_age_year = y * fleet->n_ages + a;
          fdq_["agecomp_proportion"][i_age_year] =
              fdq_["agecomp_expected"][i_age_year] / sum;
          // robust_add + robust_sum * this->agecomp_expected[i_age_year] / sum;
 
          if (fleet->fleet_observed_agecomp_data_id_m != -999) {
            fdq_["agecomp_expected"][i_age_year] =
                fdq_["agecomp_proportion"][i_age_year] * sum_obs;
          }
        }
      }
    }
  }
 
  // --- Growth-derived WAA helpers ---
  //
  // These helpers compute and cache fleet-level and population-level
  // weight-at-age values that use growth-derived age-to-length conversion rows
  // when available.

  Type MeanWeightFromAgeToLengthConversionRow(
      const std::shared_ptr<fims_popdy::GrowthDerivedObservationBase<Type>>
          &growth_observation,
      const std::shared_ptr<fims_popdy::Fleet<Type>> &fleet,
      const fims::Vector<Type> &age_to_length_conversion_row) {
    Type mean_weight = static_cast<Type>(0.0);
    for (size_t l = 0; l < fleet->n_lengths; ++l) {
      const double bin_midpoint =
          0.5 * (fleet->length_bin_edges[l] + fleet->length_bin_edges[l + 1]);
      mean_weight += age_to_length_conversion_row[l] *
                     growth_observation->EvaluateWeightAtLength(bin_midpoint);
    }
    return mean_weight;
  }

  void BuildAgeToLengthConversionDerivedRowOrThrow(
      const std::shared_ptr<fims_popdy::AgeToLengthConversionDerived<Type>>
          &age_to_length_conversion_derived,
      const std::shared_ptr<fims_popdy::Fleet<Type>> &fleet, size_t year,
      size_t age, fims::Vector<Type> &age_to_length_conversion_row) {
    if (age_to_length_conversion_derived == nullptr ||
        !age_to_length_conversion_derived->IsActive()) {
      throw std::runtime_error(
          "Growth-derived age-to-length conversion was unavailable while "
          "building fleet age-to-length probabilities.");
    }
 
    if (!age_to_length_conversion_derived->BuildAgeToLengthConversionRow(
            year, age, age_to_length_conversion_row)) {
      std::stringstream ss;
      ss << "Failed to build growth-derived age-to-length conversion row for "
            "fleet id "
         << fleet->GetId() << ", year " << year << ", age " << age << ".";
      FIMS_ERROR_LOG(ss.str());
      throw std::runtime_error(ss.str());
    }
  }

  Type MeanWeightFromAgeToLengthConversionDerived(
      const std::shared_ptr<fims_popdy::AgeToLengthConversionDerived<Type>>
          &age_to_length_conversion_derived,
      const std::shared_ptr<fims_popdy::Fleet<Type>> &fleet, size_t year,
      size_t age) {
    fims::Vector<Type> age_to_length_conversion_row;
    BuildAgeToLengthConversionDerivedRowOrThrow(
        age_to_length_conversion_derived, fleet, year, age,
        age_to_length_conversion_row);
 
    return MeanWeightFromAgeToLengthConversionRow(
        age_to_length_conversion_derived->growth_observation_, fleet,
        age_to_length_conversion_row);
  }

  void EnsureAllFleetAgeToLengthConversion() {
    for (size_t p = 0; p < this->populations.size(); ++p) {
      fims_popdy::EnsurePopulationFleetAgeToLengthConversion<Type>(
          this->populations[p]);
    }
  }

  void PreparePopulationGrowthProducts() {
    for (size_t p = 0; p < this->populations.size(); ++p) {
      std::shared_ptr<fims_popdy::Population<Type>> &population =
          this->populations[p];
 
      if (population == nullptr || population->growth == nullptr) {
        continue;
      }
 
      std::shared_ptr<fims_popdy::GrowthDerivedObservationBase<Type>>
          growth_observation = std::dynamic_pointer_cast<
              fims_popdy::GrowthDerivedObservationBase<Type>>(
              population->growth);
 
      if (growth_observation != nullptr) {
        growth_observation->PrepareGrowthProducts();
 
        // Prepared on every evaluation because nothing marks the size
        // products stale when growth parameters change.
        if (population->size_distribution_provider != nullptr) {
          population->size_distribution_provider->PrepareSizeProducts();
        }
      }
    }
  }

  Type GrowthDerivedFleetMeanWeightAA(
      const std::shared_ptr<fims_popdy::Fleet<Type>> &fleet, size_t year,
      size_t age) {
    if (fleet == nullptr) {
      throw std::runtime_error(
          "Fleet pointer was null while resolving growth-derived fleet mean "
          "weight-at-age.");
    }
 
    std::shared_ptr<fims_popdy::AgeToLengthConversionDerived<Type>>
        age_to_length_conversion_derived_model = std::dynamic_pointer_cast<
            fims_popdy::AgeToLengthConversionDerived<Type>>(
            fleet->age_to_length_conversion_model);
 
    if (age_to_length_conversion_derived_model != nullptr &&
        age_to_length_conversion_derived_model->IsActive()) {
      return MeanWeightFromAgeToLengthConversionDerived(
          age_to_length_conversion_derived_model, fleet, year, age);
    }
 
    std::stringstream ss;
    ss << "Failed to resolve growth-derived fleet mean weight-at-age for fleet "
          "id "
       << fleet->GetId() << ", year " << year << ", age " << age << ".";
    FIMS_ERROR_LOG(ss.str());
    throw std::runtime_error(ss.str());
  }

  Type BiologicalMeanWeightFromPreparedGrowthProducts(
      const std::shared_ptr<fims_popdy::Population<Type>> &population,
      const std::shared_ptr<fims_popdy::GrowthDerivedObservationBase<Type>>
          &growth_observation,
      size_t year, size_t age) {
    if (population == nullptr || growth_observation == nullptr) {
      throw std::runtime_error(
          "Population or growth-derived observation pointer was null while "
          "resolving biological mean weight-at-age.");
    }
 
    const GrowthProducts<Type> *gp =
        growth_observation->TryGetPreparedGrowthProducts();
 
    if (gp == nullptr) {
      growth_observation->PrepareGrowthProducts();
      gp = growth_observation->TryGetPreparedGrowthProducts();
    }
 
    if (gp == nullptr) {
      throw std::runtime_error(
          "Growth products were unavailable while resolving biological mean "
          "weight-at-age.");
    }
 
    if (gp->n_years == 0 || gp->n_ages == 0) {
      throw std::runtime_error(
          "Prepared growth products were empty while resolving biological mean "
          "weight-at-age.");
    }
 
    if (gp->n_sexes != 1) {
      throw std::runtime_error(
          "Biological mean weight-at-age currently requires one prepared sex.");
    }
 
    const size_t y = (std::min)(year, gp->n_years - 1);
    const size_t a = (std::min)(age, gp->n_ages - 1);
    return gp->MeanWAA(y, a, 0);
  }

  Type PopulationMeanWeightAA(
      const std::shared_ptr<fims_popdy::Population<Type>> &population,
      size_t year, size_t age) {
    if (population == nullptr) {
      throw std::runtime_error(
          "Population pointer was null while resolving population mean "
          "weight-at-age.");
    }
 
    if (population->growth == nullptr) {
      throw std::runtime_error(
          "Population growth pointer was null while resolving population mean "
          "weight-at-age.");
    }
 
    std::shared_ptr<fims_popdy::GrowthDerivedObservationBase<Type>>
        growth_observation = std::dynamic_pointer_cast<
            fims_popdy::GrowthDerivedObservationBase<Type>>(population->growth);
 
    if (growth_observation == nullptr) {
      return population->growth->evaluate(year, population->ages[age]);
    }
 
    if (population->size_distribution_provider != nullptr) {
      return population->size_distribution_provider->MeanWeightAtAge(
          (std::min)(year, population->n_years - 1), age);
    }
 
    return BiologicalMeanWeightFromPreparedGrowthProducts(
        population, growth_observation, year, age);
  }

  void evaluate_length_comp() {
    fleet_iterator fit;
    for (fit = this->fleets.begin(); fit != this->fleets.end(); ++fit) {
      std::map<std::string, fims::Vector<Type>> &fdq_ =
          this->GetFleetDerivedQuantities((*fit).second->GetId());
 
      std::shared_ptr<fims_popdy::Fleet<Type>> &fleet = (*fit).second;
 
      if (!fleet->requires_age_length_mapping || fleet->n_lengths == 0) {
        continue;
      }
 
      if (fleet->age_to_length_conversion_model == nullptr ||
          !fleet->age_to_length_conversion_model->IsActive()) {
        std::stringstream ss;
        ss << "Fleet id " << fleet->GetId()
           << "no usable age-to-length conversion path";
        FIMS_ERROR_LOG(ss.str());
        throw std::runtime_error(ss.str());
      }
 
      for (size_t y = 0; y < fleet->n_years; y++) {
        Type sum = static_cast<Type>(0.0);
        Type sum_obs = static_cast<Type>(0.0);
        // robust_add is a small value to add to expected composition
        // proportions at age to stabilize likelihood calculations
        // when the expected proportions are close to zero.
        // Type robust_add = static_cast<Type>(0.0); // 0.0001; zeroed out
        // before testing sum robust is used to calculate the total sum of
        // robust additions to ensure that proportions sum to 1. Type
        // robust_sum = static_cast<Type>(1.0);
        for (size_t a = 0; a < fleet->n_ages; a++) {
          size_t i_age_year = y * fleet->n_ages + a;
          fims::Vector<Type> age_to_length_conversion_row;
          if (!fleet->age_to_length_conversion_model
                   ->BuildAgeToLengthConversionRow(
                       y, a, age_to_length_conversion_row)) {
            std::stringstream ss;
            ss << "Failed to build age-to-length conversion row for fleet id "
               << fleet->GetId() << ", year " << y << ", age " << a << ".";
            FIMS_ERROR_LOG(ss.str());
            throw std::runtime_error(ss.str());
          }
 
          for (size_t l = 0; l < fleet->n_lengths; l++) {
            size_t i_length_year = y * fleet->n_lengths + l;
            const Type age_to_length_prob = age_to_length_conversion_row[l];
            // Use agecomp_proportion rather than agecomp_expected, which is
            // rescaled to the observed age-composition sample size and is
            // therefore 0 in years without age-composition data. The scale
            // cancels when the length composition is normalized below.
            fdq_["lengthcomp_expected"][i_length_year] +=
                fdq_["agecomp_proportion"][i_age_year] * age_to_length_prob;
 
            fdq_["catch_numbers_at_length"][i_length_year] +=
                fdq_["catch_numbers_at_age"][i_age_year] * age_to_length_prob;
 
            fdq_["index_numbers_at_length"][i_length_year] +=
                fdq_["index_numbers_at_age"][i_age_year] * age_to_length_prob;
          }
        }
 
        for (size_t l = 0; l < fleet->n_lengths; l++) {
          size_t i_length_year = y * fleet->n_lengths + l;
          sum += fdq_["lengthcomp_expected"][i_length_year];
          // robust_sum -= robust_add;
 
          if (fleet->fleet_observed_lengthcomp_data_id_m != -999) {
            if (fleet->observed_lengthcomp_data->at(i_length_year) !=
                fleet->observed_lengthcomp_data->na_value) {
              sum_obs += fleet->observed_lengthcomp_data->at(i_length_year);
            }
          }
        }
        for (size_t l = 0; l < fleet->n_lengths; l++) {
          size_t i_length_year = y * fleet->n_lengths + l;
          fdq_["lengthcomp_proportion"][i_length_year] =
              fdq_["lengthcomp_expected"][i_length_year] / sum;
          // robust_add + robust_sum *
          // this->lengthcomp_expected[i_length_year] / sum;
          if (fleet->fleet_observed_lengthcomp_data_id_m != -999) {
            fdq_["lengthcomp_expected"][i_length_year] =
                fdq_["lengthcomp_proportion"][i_length_year] * sum_obs;
          }
        }
      }
    }
  }

  void evaluate_ageing_error() {
    for (fleet_iterator fit = this->fleets.begin(); fit != this->fleets.end();
         ++fit) {
      std::shared_ptr<fims_popdy::Fleet<Type>> &fleet = (*fit).second;
      if (fleet->ageing_error.size() == 0) {
        continue;
      }
      std::map<std::string, fims::Vector<Type>> &fdq_ =
          this->GetFleetDerivedQuantities(fleet->GetId());
 
      fims::Vector<Type> proportion_at_true_age(fleet->n_ages);
      fims::Vector<Type> expected_at_true_age(fleet->n_ages);
      for (size_t y = 0; y < fleet->n_years; y++) {
        for (size_t a = 0; a < fleet->n_ages; a++) {
          size_t i_age_year = y * fleet->n_ages + a;
          proportion_at_true_age[a] = fdq_["agecomp_proportion"][i_age_year];
          expected_at_true_age[a] = fdq_["agecomp_expected"][i_age_year];
        }
        fims::Vector<Type> proportion_at_observed_age =
            fims_math::apply_ageing_error(proportion_at_true_age,
                                          fleet->ageing_error, y);
        fims::Vector<Type> expected_at_observed_age =
            fims_math::apply_ageing_error(expected_at_true_age,
                                          fleet->ageing_error, y);
        for (size_t a = 0; a < fleet->n_ages; a++) {
          size_t i_age_year = y * fleet->n_ages + a;
          fdq_["agecomp_proportion"][i_age_year] =
              proportion_at_observed_age[a];
          fdq_["agecomp_expected"][i_age_year] = expected_at_observed_age[a];
        }
      }
    }
  }

  void evaluate_index() {
    fleet_iterator fit;
    for (fit = this->fleets.begin(); fit != this->fleets.end(); ++fit) {
      std::map<std::string, fims::Vector<Type>> &fdq_ =
          this->GetFleetDerivedQuantities((*fit).second->GetId());
      std::shared_ptr<fims_popdy::Fleet<Type>> &fleet = (*fit).second;
 
      for (size_t i = 0; i < fdq_["index_numbers"].size(); i++) {
        if (fleet->observed_index_units == "number") {
          fdq_["index_expected"][i] = fdq_["index_numbers"][i];
        } else {
          fdq_["index_expected"][i] = fdq_["index_weight"][i];
        }
        fdq_["log_index_expected"][i] = log(fdq_["index_expected"][i]);
      }
    }
  }

  void evaluate_catch() {
    fleet_iterator fit;
    for (fit = this->fleets.begin(); fit != this->fleets.end(); ++fit) {
      std::map<std::string, fims::Vector<Type>> &fdq_ =
          this->GetFleetDerivedQuantities((*fit).second->GetId());
      std::shared_ptr<fims_popdy::Fleet<Type>> &fleet = (*fit).second;
 
      for (size_t i = 0; i < fdq_["catch_weight"].size(); i++) {
        if (fleet->observed_catch_units == "number") {
          fdq_["catch_expected"][i] = fdq_["catch_numbers"][i];
        } else {
          fdq_["catch_expected"][i] = fdq_["catch_weight"][i];
        }
        fdq_["log_catch_expected"][i] = log(fdq_["catch_expected"][i]);
      }
    }
  }
 
  virtual void Evaluate() {
    /*
               Sets derived vectors to zero
               Performs parameters transformations
               Sets recruitment deviations to mean 0.
     */
    Prepare();
    PreparePopulationGrowthProducts();
    EnsureAllFleetAgeToLengthConversion();
    /*
     start at year=0, age=0;
     here year 0 is the estimated initial population structure and age 0 are
     recruits loops start at zero with if statements inside to specify unique
     code for initial structure and recruitment 0 loops. Could also have started
     loops at 1 with initial structure and recruitment setup outside the loops.
 
     year loop is extended to <= n_years because SB is calculated as the start
     of the year value and by extending one extra year we get estimates of the
     population structure at the end of the final year. An alternative approach
     would be to keep initial numbers at age in it's own vector and each year to
     include the population structure at the end of the year. This is likely a
     null point given that we are planning to modify to an event/stanza based
     structure in later milestones which will eliminate this confusion by
     explicitly referencing the exact date (or period of averaging) at which any
     calculation or output is being made.
     */
    for (size_t p = 0; p < this->populations.size(); p++) {
      std::shared_ptr<fims_popdy::Population<Type>> &population =
          this->populations[p];
      std::map<std::string, fims::Vector<Type>> &pdq_ =
          this->GetPopulationDerivedQuantities(population->GetId());
      // CAAPopulationProxy<Type>& population = this->populations_proxies[p];
 
      for (size_t y = 0; y <= population->n_years; y++) {
        for (size_t a = 0; a < population->n_ages; a++) {
          /*
           index naming defines the dimensional folding structure
           i.e. i_age_year is referencing folding over years and ages.
           */
          size_t i_age_year = y * population->n_ages + a;
          /*
           Mortality rates are not estimated in the final year which is
           used to show expected population structure at the end of the model
           period. This is because biomass in year i represents biomass at the
           start of the year. Should we add complexity to track more values such
           as start, mid, and end biomass in all years where, start biomass=end
           biomass of the previous year? Referenced above, this is probably not
           worth exploring as later milestone changes will eliminate this
           confusion.
           */
          if (y < population->n_years) {
            /*
             First thing we need is total mortality aggregated across all fleets
             to inform the subsequent catch and change in numbers at age
             calculations. This is only calculated for years < n_years as these
             are the model estimated years with data. The year loop extends to
             y=n_years so that population numbers at age and SB can be
             calculated at the end of the last year of the model
             */
            CalculateMortality(population, i_age_year, y, a);
          }
          CalculateMaturityAA(population, i_age_year, y, a);
          /* if statements needed because some quantities are only needed
          for the first year and/or age, so these steps are included here.
           */
          if (y == 0) {
            // Initial numbers at age is a user input or estimated parameter
            // vector.
            CalculateInitialNumbersAA(population, i_age_year, a);
 
            if (a == 0) {
              /*
             Expected recruitment in year 0 is numbers at age 0 in year 0.
             */
              pdq_["expected_recruitment"][y] =
                  pdq_["numbers_at_age"][i_age_year];
              pdq_["unfished_numbers_at_age"][i_age_year] =
                  fims_math::exp(population->recruitment->log_rzero[0]);
            } else {
              CalculateUnfishedInitialNumbersAA(population, i_age_year, a);
            }
 
          } else {
            if (a == 0) {
              // Set the nrecruits for age a=0 year y (use pointers instead of
              // functional returns)
              CalculateRecruitment(population, i_age_year, y, y);
              pdq_["unfished_numbers_at_age"][i_age_year] =
                  fims_math::exp(population->recruitment->log_rzero[0]);
            } else {
              size_t i_agem1_yearm1 = (y - 1) * population->n_ages + (a - 1);
              CalculateNumbersAA(population, i_age_year, i_agem1_yearm1, a);
              CalculateUnfishedNumbersAA(population, i_age_year, i_agem1_yearm1,
                                         a);
            }
          }
 
          /*
            Fished and unfished biomass vectors are summing biomass at
            age across ages.
          */
 
          CalculateBiomass(population, i_age_year, y, a);
 
          CalculateUnfishedBiomass(population, i_age_year, y, a);
 
          /*
            Fished and unfished spawning biomass vectors are summing biomass at
            age across ages to allow calculation of recruitment in the next
            year.
          */
 
          CalculateSpawningBiomass(population, i_age_year, y, a);
 
          CalculateUnfishedSpawningBiomass(population, i_age_year, y, a);
 
          /*
          Here composition, total catch, and index values are calculated for all
          years with reference data. They are not calculated for y=n_years as
          there is this is just to get final population structure at the end of
          the terminal year.
           */
          if (y < population->n_years) {
            CalculateCatchNumbersAA(population, i_age_year, y, a);
            CalculateCatchWeightAA(population, y, a);
            CalculateCatch(population, y, a);
 
            CalculateIndexNumbersAA(population, i_age_year, y, a);
            CalculateIndexWeightAA(population, y, a);
            CalculateIndex(population, i_age_year, y, a);
          }
        }
        /* Calculate spawning biomass depletion ratio */
        CalculateSpawningBiomassRatio(population, y);
      }
    }
    evaluate_age_comp();
    evaluate_length_comp();
    evaluate_ageing_error();
    evaluate_index();
    evaluate_catch();
  }

  virtual void Report() {
    int n_fleets = this->fleets.size();
    int n_pops = this->populations.size();
#ifdef TMB_MODEL
    if (this->do_reporting == true) {
      EnsureAllFleetAgeToLengthConversion();
      report_vectors.clear();
      // std::shared_ptr<UncertaintyReportInfoMap>
      // population_uncertainty_report_info_map =
      //     this->GetPopulationUncertaintyReportInfoMap();
 
      // std::shared_ptr<UncertaintyReportInfoMap>
      // fleet_uncertainty_report_info_map =
      //     this->GetFleetUncertaintyReportInfoMap();
 
      // initialize population vectors
      vector<vector<Type>> biomass_p(n_pops);
      vector<vector<Type>> expected_recruitment_p(n_pops);
      vector<vector<Type>> mortality_F_p(n_pops);
      vector<vector<Type>> mortality_M_p(n_pops);
      vector<vector<Type>> mortality_Z_p(n_pops);
      vector<vector<Type>> numbers_at_age_p(n_pops);
      vector<vector<Type>> proportion_mature_at_age_p(n_pops);
      vector<vector<Type>> spawning_biomass_p(n_pops);
      vector<vector<Type>> sum_selectivity_p(n_pops);
      vector<vector<Type>> total_catch_numbers_p(n_pops);
      vector<vector<Type>> total_catch_weight_p(n_pops);
      vector<vector<Type>> unfished_biomass_p(n_pops);
      vector<vector<Type>> unfished_numbers_at_age_p(n_pops);
      vector<vector<Type>> unfished_spawning_biomass_p(n_pops);
      vector<vector<Type>> spawning_biomass_ratio_p(n_pops);
      vector<vector<Type>> growth_mean_LAA_p(n_pops);
      vector<vector<Type>> growth_sd_LAA_p(n_pops);
      vector<vector<Type>> growth_mean_WAA_p(n_pops);
 
      // initialize fleet vectors
      vector<vector<Type>> agecomp_expected_f(n_fleets);
      vector<vector<Type>> agecomp_proportion_f(n_fleets);
      vector<vector<Type>> age_to_length_conversion_f(n_fleets);
      vector<vector<Type>> age_to_length_conversion_derived_f(n_fleets);
      vector<vector<Type>> age_to_length_conversion_derived_used_f(n_fleets);
      vector<vector<Type>> growth_derived_mean_WAA_f(n_fleets);
      vector<vector<Type>> index_expected_f(n_fleets);
      vector<vector<Type>> index_numbers_f(n_fleets);
      vector<vector<Type>> index_numbers_at_age_f(n_fleets);
      vector<vector<Type>> index_numbers_at_length_f(n_fleets);
      vector<vector<Type>> index_weight_f(n_fleets);
      vector<vector<Type>> index_weight_at_age_f(n_fleets);
      vector<vector<Type>> catch_expected_f(n_fleets);
      vector<vector<Type>> catch_numbers_f(n_fleets);
      vector<vector<Type>> catch_numbers_at_age_f(n_fleets);
      vector<vector<Type>> catch_numbers_at_length_f(n_fleets);
      vector<vector<Type>> catch_weight_f(n_fleets);
      vector<vector<Type>> catch_weight_at_age_f(n_fleets);
      vector<vector<Type>> lengthcomp_expected_f(n_fleets);
      vector<vector<Type>> lengthcomp_proportion_f(n_fleets);
      vector<vector<Type>> log_index_expected_f(n_fleets);
      vector<vector<Type>> log_catch_expected_f(n_fleets);
 
      // initiate population index for structuring report out objects
      int pop_idx = 0;
      for (size_t p = 0; p < this->populations.size(); p++) {
        std::map<std::string, fims::Vector<Type>> &derived_quantities =
            this->GetPopulationDerivedQuantities(this->populations[p]->GetId());
        biomass_p(pop_idx) = derived_quantities["biomass"].to_tmb();
        expected_recruitment_p(pop_idx) =
            derived_quantities["expected_recruitment"].to_tmb();
        mortality_F_p(pop_idx) = derived_quantities["mortality_F"].to_tmb();
        mortality_M_p(pop_idx) = derived_quantities["mortality_M"].to_tmb();
        mortality_Z_p(pop_idx) = derived_quantities["mortality_Z"].to_tmb();
        numbers_at_age_p(pop_idx) =
            derived_quantities["numbers_at_age"].to_tmb();
        proportion_mature_at_age_p(pop_idx) =
            derived_quantities["proportion_mature_at_age"].to_tmb();
        spawning_biomass_p(pop_idx) =
            derived_quantities["spawning_biomass"].to_tmb();
        sum_selectivity_p(pop_idx) =
            derived_quantities["sum_selectivity"].to_tmb();
        total_catch_numbers_p(pop_idx) =
            derived_quantities["total_catch_numbers"].to_tmb();
        total_catch_weight_p(pop_idx) =
            derived_quantities["total_catch_weight"].to_tmb();
        unfished_biomass_p(pop_idx) =
            derived_quantities["unfished_biomass"].to_tmb();
        unfished_numbers_at_age_p(pop_idx) =
            derived_quantities["unfished_numbers_at_age"].to_tmb();
        unfished_spawning_biomass_p(pop_idx) =
            derived_quantities["unfished_spawning_biomass"].to_tmb();
        spawning_biomass_ratio_p(pop_idx) =
            this->populations[pop_idx]->spawning_biomass_ratio.to_tmb();
 
        if (std::shared_ptr<fims_popdy::GrowthDerivedObservationBase<Type>>
                growth_observation = std::dynamic_pointer_cast<
                    fims_popdy::GrowthDerivedObservationBase<Type>>(
                    this->populations[p]->growth)) {
          const auto *gp = growth_observation->TryGetPreparedGrowthProducts();
          if (gp == nullptr) {
            throw std::runtime_error(
                "Growth products were not prepared before report generation.");
          }
 
          const std::size_t n = gp->Size();
          vector<Type> mean_laa(n);
          vector<Type> sd_laa(n);
          vector<Type> mean_waa(n);
          for (std::size_t i = 0; i < n; ++i) {
            mean_laa(i) = gp->mean_LAA[i];
            sd_laa(i) = gp->sd_LAA[i];
            const std::size_t year = i / this->populations[p]->n_ages;
            const std::size_t age = i % this->populations[p]->n_ages;
            mean_waa(i) =
                PopulationMeanWeightAA(this->populations[p], year, age);
          }
          growth_mean_LAA_p(pop_idx) = mean_laa;
          growth_sd_LAA_p(pop_idx) = sd_laa;
          growth_mean_WAA_p(pop_idx) = mean_waa;
        }
 
        pop_idx += 1;
      }
 
      // initiate fleet index for structuring report out objects
      int fleet_idx = 0;
      fleet_iterator fit;
      for (fit = this->fleets.begin(); fit != this->fleets.end(); ++fit) {
        std::shared_ptr<fims_popdy::Fleet<Type>> &fleet = (*fit).second;
        std::map<std::string, fims::Vector<Type>> &derived_quantities =
            this->GetFleetDerivedQuantities(fleet->GetId());
 
        agecomp_expected_f(fleet_idx) =
            derived_quantities["agecomp_expected"].to_tmb();
        agecomp_proportion_f(fleet_idx) =
            derived_quantities["agecomp_proportion"].to_tmb();
        fims::Vector<Type> age_to_length_conversion_derived;
        fims::Vector<Type> age_to_length_conversion_derived_used;
        fims::Vector<Type> growth_derived_mean_WAA;
        bool using_age_to_length_conversion_derived = false;
        if (isDouble<Type>::value) {
          age_to_length_conversion_derived_used.resize(1);
          age_to_length_conversion_derived_used[0] = static_cast<Type>(0.0);
 
          std::shared_ptr<fims_popdy::AgeToLengthConversionDerived<Type>>
              age_to_length_conversion_derived_model =
                  std::dynamic_pointer_cast<
                      fims_popdy::AgeToLengthConversionDerived<Type>>(
                      fleet->age_to_length_conversion_model);
 
          if (age_to_length_conversion_derived_model != nullptr &&
              age_to_length_conversion_derived_model->IsActive()) {
            using_age_to_length_conversion_derived = true;
            age_to_length_conversion_derived_used[0] = static_cast<Type>(1.0);
            growth_derived_mean_WAA.resize(fleet->n_years * fleet->n_ages);
 
            for (size_t y = 0; y < fleet->n_years; ++y) {
              for (size_t a = 0; a < fleet->n_ages; ++a) {
                const size_t i_age_year = y * fleet->n_ages + a;
 
                if (this->report_age_to_length_conversion_derived_tensor) {
                  if (age_to_length_conversion_derived.size() == 0) {
                    age_to_length_conversion_derived.resize(
                        fleet->n_years * fleet->n_ages * fleet->n_lengths);
                  }
 
                  fims::Vector<Type> age_to_length_conversion_row;
                  BuildAgeToLengthConversionDerivedRowOrThrow(
                      age_to_length_conversion_derived_model, fleet, y, a,
                      age_to_length_conversion_row);
 
                  growth_derived_mean_WAA[i_age_year] =
                      MeanWeightFromAgeToLengthConversionRow(
                          age_to_length_conversion_derived_model
                              ->growth_observation_,
                          fleet, age_to_length_conversion_row);
 
                  for (size_t l = 0; l < fleet->n_lengths; ++l) {
                    const size_t i_length_age_year =
                        y * (fleet->n_ages * fleet->n_lengths) +
                        a * fleet->n_lengths + l;
                    age_to_length_conversion_derived[i_length_age_year] =
                        age_to_length_conversion_row[l];
                  }
                } else {
                  growth_derived_mean_WAA[i_age_year] =
                      GrowthDerivedFleetMeanWeightAA(fleet, y, a);
                }
              }
            }
          }
        }
 
        age_to_length_conversion_derived_f(fleet_idx) =
            age_to_length_conversion_derived.to_tmb();
        age_to_length_conversion_derived_used_f(fleet_idx) =
            age_to_length_conversion_derived_used.to_tmb();
        growth_derived_mean_WAA_f(fleet_idx) = growth_derived_mean_WAA.to_tmb();
        if (!using_age_to_length_conversion_derived) {
          age_to_length_conversion_f(fleet_idx) =
              fleet->age_to_length_conversion.to_tmb();
        }
        index_expected_f(fleet_idx) =
            derived_quantities["index_expected"].to_tmb();
        index_numbers_f(fleet_idx) =
            derived_quantities["index_numbers"].to_tmb();
        index_numbers_at_age_f(fleet_idx) =
            derived_quantities["index_numbers_at_age"].to_tmb();
        index_numbers_at_length_f(fleet_idx) =
            derived_quantities["index_numbers_at_length"].to_tmb();
        index_weight_f(fleet_idx) = derived_quantities["index_weight"].to_tmb();
        index_weight_at_age_f(fleet_idx) =
            derived_quantities["index_weight_at_age"].to_tmb();
        catch_expected_f(fleet_idx) =
            derived_quantities["catch_expected"].to_tmb();
        catch_numbers_f(fleet_idx) =
            derived_quantities["catch_numbers"].to_tmb();
        catch_numbers_at_age_f(fleet_idx) =
            derived_quantities["catch_numbers_at_age"].to_tmb();
        catch_numbers_at_length_f(fleet_idx) =
            derived_quantities["catch_numbers_at_length"].to_tmb();
        catch_weight_f(fleet_idx) = derived_quantities["catch_weight"].to_tmb();
        catch_weight_at_age_f(fleet_idx) =
            derived_quantities["catch_weight_at_age"].to_tmb();
        // length_comp_expected_f(fleet_idx) =
        // derived_quantities["length_comp_expected"];
        // length_comp_proportion_f(fleet_idx) =
        // derived_quantities["length_comp_proportion"];
        lengthcomp_expected_f(fleet_idx) =
            derived_quantities["lengthcomp_expected"].to_tmb();
        lengthcomp_proportion_f(fleet_idx) =
            derived_quantities["lengthcomp_proportion"].to_tmb();
        log_index_expected_f(fleet_idx) =
            derived_quantities["log_index_expected"].to_tmb();
        log_catch_expected_f(fleet_idx) =
            derived_quantities["log_catch_expected"].to_tmb();
        fleet_idx += 1;
      }
 
      vector<Type> biomass = ADREPORTvector(biomass_p);
      vector<Type> expected_recruitment =
          ADREPORTvector(expected_recruitment_p);
      vector<Type> mortality_F = ADREPORTvector(mortality_F_p);
      vector<Type> mortality_M = ADREPORTvector(mortality_M_p);
      vector<Type> mortality_Z = ADREPORTvector(mortality_Z_p);
      vector<Type> numbers_at_age = ADREPORTvector(numbers_at_age_p);
      vector<Type> proportion_mature_at_age =
          ADREPORTvector(proportion_mature_at_age_p);
      vector<Type> spawning_biomass = ADREPORTvector(spawning_biomass_p);
      vector<Type> sum_selectivity = ADREPORTvector(sum_selectivity_p);
      vector<Type> total_catch_numbers = ADREPORTvector(total_catch_numbers_p);
      vector<Type> total_catch_weight = ADREPORTvector(total_catch_weight_p);
      vector<Type> unfished_biomass = ADREPORTvector(unfished_biomass_p);
      vector<Type> unfished_numbers_at_age =
          ADREPORTvector(unfished_numbers_at_age_p);
      vector<Type> unfished_spawning_biomass =
          ADREPORTvector(unfished_spawning_biomass_p);
      vector<Type> spawning_biomass_ratio =
          ADREPORTvector(spawning_biomass_ratio_p);
 
      vector<Type> agecomp_expected = ADREPORTvector(agecomp_expected_f);
      vector<Type> agecomp_proportion = ADREPORTvector(agecomp_proportion_f);
      vector<Type> index_expected = ADREPORTvector(index_expected_f);
      vector<Type> index_numbers = ADREPORTvector(index_numbers_f);
      vector<Type> index_numbers_at_age =
          ADREPORTvector(index_numbers_at_age_f);
      vector<Type> index_numbers_at_length =
          ADREPORTvector(index_numbers_at_length_f);
      vector<Type> index_weight = ADREPORTvector(index_weight_f);
      vector<Type> index_weight_at_age = ADREPORTvector(index_weight_at_age_f);
      vector<Type> catch_expected = ADREPORTvector(catch_expected_f);
      vector<Type> catch_numbers = ADREPORTvector(catch_numbers_f);
      vector<Type> catch_numbers_at_age =
          ADREPORTvector(catch_numbers_at_age_f);
      vector<Type> catch_numbers_at_length =
          ADREPORTvector(catch_numbers_at_length_f);
      vector<Type> catch_weight = ADREPORTvector(catch_weight_f);
      vector<Type> catch_weight_at_age = ADREPORTvector(catch_weight_at_age_f);
      // vector<Type> length_comp_expected =
      // ADREPORTvector(length_comp_expected_f); vector<Type>
      // length_comp_proportion = ADREPORTvector(length_comp_proportion_f);
      vector<Type> lengthcomp_expected = ADREPORTvector(lengthcomp_expected_f);
      vector<Type> lengthcomp_proportion =
          ADREPORTvector(lengthcomp_proportion_f);
      vector<Type> log_index_expected = ADREPORTvector(log_index_expected_f);
      vector<Type> log_catch_expected = ADREPORTvector(log_catch_expected_f);
      // populations
      // report
      FIMS_REPORT_F_("biomass", biomass_p, this->of);
      FIMS_REPORT_F_("expected_recruitment", expected_recruitment_p, this->of);
      FIMS_REPORT_F_("mortality_F", mortality_F_p, this->of);
      FIMS_REPORT_F_("mortality_M", mortality_M_p, this->of);
      FIMS_REPORT_F_("mortality_Z", mortality_Z_p, this->of);
      FIMS_REPORT_F_("numbers_at_age", numbers_at_age_p, this->of);
      FIMS_REPORT_F_("proportion_mature_at_age", proportion_mature_at_age_p,
                     this->of);
      FIMS_REPORT_F_("spawning_biomass", spawning_biomass_p, this->of);
      FIMS_REPORT_F_("sum_selectivity", sum_selectivity_p, this->of);
      FIMS_REPORT_F_("total_catch_numbers", total_catch_numbers_p, this->of);
      FIMS_REPORT_F_("total_catch_weight", total_catch_weight_p, this->of);
      FIMS_REPORT_F_("unfished_biomass", unfished_biomass_p, this->of);
      FIMS_REPORT_F_("unfished_numbers_at_age", unfished_numbers_at_age_p,
                     this->of);
      FIMS_REPORT_F_("unfished_spawning_biomass", unfished_spawning_biomass_p,
                     this->of);
      FIMS_REPORT_F_("spawning_biomass_ratio", spawning_biomass_ratio_p,
                     this->of);
      FIMS_REPORT_F_("growth_mean_LAA", growth_mean_LAA_p, this->of);
      FIMS_REPORT_F_("growth_sd_LAA", growth_sd_LAA_p, this->of);
      FIMS_REPORT_F_("growth_mean_WAA", growth_mean_WAA_p, this->of);
 
      // adreport
      ADREPORT_F(biomass, this->of);
      ADREPORT_F(expected_recruitment, this->of);
      ADREPORT_F(mortality_F, this->of);
      ADREPORT_F(mortality_M, this->of);
      ADREPORT_F(mortality_Z, this->of);
      ADREPORT_F(numbers_at_age, this->of);
      ADREPORT_F(proportion_mature_at_age, this->of);
      ADREPORT_F(spawning_biomass, this->of);
      ADREPORT_F(sum_selectivity, this->of);
      ADREPORT_F(total_catch_numbers, this->of);
      ADREPORT_F(total_catch_weight, this->of);
      ADREPORT_F(unfished_biomass, this->of);
      ADREPORT_F(unfished_numbers_at_age, this->of);
      ADREPORT_F(unfished_spawning_biomass, this->of);
      ADREPORT_F(spawning_biomass_ratio, this->of);
 
      // fleets
      // report
      FIMS_REPORT_F_("agecomp_expected", agecomp_expected_f, this->of);
      FIMS_REPORT_F_("agecomp_proportion", agecomp_proportion_f, this->of);
      FIMS_REPORT_F_("age_to_length_conversion", age_to_length_conversion_f,
                     this->of);
      FIMS_REPORT_F_("age_to_length_conversion_derived",
                     age_to_length_conversion_derived_f, this->of);
      FIMS_REPORT_F_("age_to_length_conversion_derived_used",
                     age_to_length_conversion_derived_used_f, this->of);
      FIMS_REPORT_F_("growth_derived_mean_WAA", growth_derived_mean_WAA_f,
                     this->of);
      FIMS_REPORT_F_("index_expected", index_expected_f, this->of);
      FIMS_REPORT_F_("index_numbers", index_numbers_f, this->of);
      FIMS_REPORT_F_("index_numbers_at_age", index_numbers_at_age_f, this->of);
      FIMS_REPORT_F_("index_numbers_at_length", index_numbers_at_length_f,
                     this->of);
      FIMS_REPORT_F_("index_weight", index_weight_f, this->of);
      FIMS_REPORT_F_("index_weight_at_age", index_weight_at_age_f, this->of);
      FIMS_REPORT_F_("catch_expected", catch_expected_f, this->of);
      FIMS_REPORT_F_("catch_numbers", catch_numbers_f, this->of);
      FIMS_REPORT_F_("catch_numbers_at_age", catch_numbers_at_age_f, this->of);
      FIMS_REPORT_F_("catch_numbers_at_length", catch_numbers_at_length_f,
                     this->of);
      FIMS_REPORT_F_("catch_weight", catch_weight_f, this->of);
      FIMS_REPORT_F_("catch_weight_at_age", catch_weight_at_age_f, this->of);
      FIMS_REPORT_F_("lengthcomp_expected", lengthcomp_expected_f, this->of);
      FIMS_REPORT_F_("lengthcomp_proportion", lengthcomp_proportion_f,
                     this->of);
      FIMS_REPORT_F_("log_index_expected", log_index_expected_f, this->of);
      FIMS_REPORT_F_("log_catch_expected", log_catch_expected_f, this->of);
      // adreport
      ADREPORT_F(agecomp_expected, this->of);
      ADREPORT_F(agecomp_proportion, this->of);
      ADREPORT_F(index_expected, this->of);
      ADREPORT_F(index_numbers, this->of);
      ADREPORT_F(index_numbers_at_age, this->of);
      ADREPORT_F(index_numbers_at_length, this->of);
      ADREPORT_F(index_weight, this->of);
      ADREPORT_F(index_weight_at_age, this->of);
      ADREPORT_F(catch_expected, this->of);
      ADREPORT_F(catch_numbers, this->of);
      ADREPORT_F(catch_numbers_at_age, this->of);
      ADREPORT_F(catch_numbers_at_length, this->of);
      ADREPORT_F(catch_weight, this->of);
      ADREPORT_F(catch_weight_at_age, this->of);
      ADREPORT_F(lengthcomp_expected, this->of);
      ADREPORT_F(lengthcomp_proportion, this->of);
      ADREPORT_F(log_index_expected, this->of);
      ADREPORT_F(log_catch_expected, this->of);
      std::stringstream var_name;
      typename std::map<std::string, fims::Vector<fims::Vector<Type>>>::iterator
          rvit;
      for (rvit = report_vectors.begin(); rvit != report_vectors.end();
           ++rvit) {
        auto &x = rvit->second;
 
        int outer_dim = x.size();
        int dim = 0;
        for (int i = 0; i < outer_dim; i++) {
          dim += x[i].size();
        }
        vector<Type> res(dim);
        int idx = 0;
        for (int i = 0; i < outer_dim; i++) {
          int inner_dim = x[i].size();
          for (int j = 0; j < inner_dim; j++) {
            res(idx) = x[i][j];
            idx += 1;
          }
        }
        this->of->reportvector.push(res, rvit->first.c_str());
      }
    }
#endif
  }
};
 
}  // namespace fims_popdy
 
#endif