growth_derived_size_provider.hpp

Go to the documentation of this file.

#ifndef FIMS_POPULATION_DYNAMICS_GROWTH_DERIVED_SIZE_PROVIDER_HPP
#define FIMS_POPULATION_DYNAMICS_GROWTH_DERIVED_SIZE_PROVIDER_HPP
 
#include <cstddef>
#include <memory>
#include <stdexcept>
 
#include "../growth/growth_model_adapter.hpp"
#include "size_distribution_provider_base.hpp"
#include "size_products.hpp"
#include "../../common/def.hpp"
#include "../../common/fims_math.hpp"
 
namespace fims_popdy {

template <typename Type>
class GrowthDerivedSizeProvider : public SizeDistributionProviderBase<Type> {
 public:
  GrowthDerivedSizeProvider() = default;

  explicit GrowthDerivedSizeProvider(
      std::shared_ptr<GrowthDerivedObservationBase<Type>> growth)
      : growth_observation_(growth) {}
 
  const SizeGrid* TryGetSizeGrid() const override {
    return population_size_grid_;
  }
 
  void SetPopulationSizeGrid(const SizeGrid* population_size_grid) override {
    population_size_grid_ = population_size_grid;
    size_products_prepared_ = false;
    plus_group_warning_emitted_ = false;
  }
 
  void SetPopulationDimensions(std::size_t n_years,
                               std::size_t n_ages) override {
    n_years_ = n_years;
    n_ages_ = n_ages;
    size_products_prepared_ = false;
    plus_group_warning_emitted_ = false;
  }
 
  void PrepareSizeProducts() override {
    if (!population_size_grid_) {
      throw std::runtime_error(
          "GrowthDerivedSizeProvider requires a population biological size "
          "grid");
    }
 
    if (!population_size_grid_->IsConsistent()) {
      throw std::runtime_error(
          "GrowthDerivedSizeProvider requires a consistent population "
          "biological size grid");
    }
 
    if (population_size_grid_->n_bins == 0) {
      throw std::runtime_error(
          "GrowthDerivedSizeProvider requires a population biological size grid"
          " with at least one bin");
    }
 
    if (n_years_ == 0 || n_ages_ == 0) {
      throw std::runtime_error(
          "GrowthDerivedSizeProvider requires nonzero population dimensions");
    }
 
    const GrowthProducts<Type>& growth_products = PreparedGrowthProducts();
 
    if (growth_products.n_years != n_years_ ||
        growth_products.n_ages != n_ages_) {
      throw std::runtime_error(
          "GrowthDerivedSizeProvider growth products do not match "
          "configured population dimensions");
    }
 
    if (growth_products.n_sexes != 1) {
      throw std::runtime_error(
          "GrowthDerivedSizeProvider currently requires pooled growth "
          "products with n_sexes == 1");
    }
 
    size_products_.ResizeProbSizeOnly(n_years_, n_ages_,
                                      population_size_grid_->n_bins);
 
    const Type plus_group_warning_threshold = static_cast<Type>(0.01);
    bool plus_group_warning_needed = false;
    std::size_t warning_year_index = 0;
    std::size_t warning_age_index = 0;
 
    for (std::size_t year_index = 0; year_index < n_years_; ++year_index) {
      for (std::size_t age_index = 0; age_index < n_ages_; ++age_index) {
        fims::Vector<Type> prob_row(population_size_grid_->n_bins);
        Type row_sum = static_cast<Type>(0.0);
        const Type minimum_bin_prob = static_cast<Type>(1e-12);
 
        for (std::size_t size_bin_index = 0;
             size_bin_index < population_size_grid_->n_bins; ++size_bin_index) {
          Type bin_prob = PopulationSizeBinProb(growth_products, year_index,
                                                age_index, size_bin_index) +
                          minimum_bin_prob;
 
          prob_row[size_bin_index] = bin_prob;
          row_sum += bin_prob;
        }
 
        for (std::size_t size_bin_index = 0;
             size_bin_index < population_size_grid_->n_bins; ++size_bin_index) {
          size_products_.ProbSize(year_index, age_index, size_bin_index) =
              prob_row[size_bin_index] / row_sum;
        }
 
        const Type plus_group_prob = size_products_.ProbSize(
            year_index, age_index, population_size_grid_->n_bins - 1);
 
        if (!plus_group_warning_needed &&
            plus_group_prob > plus_group_warning_threshold) {
          plus_group_warning_needed = true;
          warning_year_index = year_index;
          warning_age_index = age_index;
        }
      }
    }
 
    if (plus_group_warning_needed && !plus_group_warning_emitted_) {
      FIMS_WARNING_LOG(
          "Growth-derived size preparation placed more than 0.01 "
          "probability in the terminal biological plus-group bin for at "
          "least one year-age row. First occurrence was year index " +
          fims::to_string(warning_year_index) + ", age index " +
          fims::to_string(warning_age_index) +
          ". The population biological size grid may need to be widened.");
      plus_group_warning_emitted_ = true;
    }
 
    // Weight depends only on length, so it is evaluated once per bin. The
    // first and last bins are weighted at their centers even though they hold
    // the lower and upper tails, as in SS3.
    fims::Vector<Type> weight_at_bin(population_size_grid_->n_bins);
    for (std::size_t size_bin_index = 0;
         size_bin_index < population_size_grid_->n_bins; ++size_bin_index) {
      weight_at_bin[size_bin_index] =
          growth_observation_->EvaluateWeightAtLength(static_cast<Type>(
              population_size_grid_->centers[size_bin_index]));
    }
 
    mean_weight_at_age_.resize(n_years_ * n_ages_);
    for (std::size_t year_index = 0; year_index < n_years_; ++year_index) {
      for (std::size_t age_index = 0; age_index < n_ages_; ++age_index) {
        Type mean_weight = static_cast<Type>(0.0);
        for (std::size_t size_bin_index = 0;
             size_bin_index < population_size_grid_->n_bins; ++size_bin_index) {
          mean_weight +=
              size_products_.ProbSize(year_index, age_index, size_bin_index) *
              weight_at_bin[size_bin_index];
        }
        mean_weight_at_age_[year_index * n_ages_ + age_index] = mean_weight;
      }
    }
 
    size_products_prepared_ = true;
  }
 
  void InvalidatePreparedSizeProducts() override {
    size_products_prepared_ = false;
  }
 
  const Type& MeanLAA(std::size_t year_index,
                      std::size_t age_index) const override {
    return PreparedGrowthProducts().MeanLAA(year_index, age_index, 0);
  }
 
  const Type& SdLAA(std::size_t year_index,
                    std::size_t age_index) const override {
    return PreparedGrowthProducts().SdLAA(year_index, age_index, 0);
  }
 
  const Type& ProbSize(std::size_t year_index, std::size_t age_index,
                       std::size_t size_bin_index) const override {
    if (!size_products_prepared_) {
      throw std::runtime_error(
          "GrowthDerivedSizeProvider requires prepared size products");
    }
 
    return size_products_.ProbSize(year_index, age_index, size_bin_index);
  }

  const Type& MeanWeightAtAge(std::size_t year_index,
                              std::size_t age_index) const override {
    if (!size_products_prepared_) {
      throw std::runtime_error(
          "GrowthDerivedSizeProvider requires prepared size products");
    }
 
    if (year_index >= n_years_ || age_index >= n_ages_) {
      throw std::out_of_range(
          "GrowthDerivedSizeProvider mean weight-at-age index out of range");
    }
 
    return mean_weight_at_age_[year_index * n_ages_ + age_index];
  }

  void SetGrowth(std::shared_ptr<GrowthDerivedObservationBase<Type>> growth) {
    growth_observation_ = growth;
    size_products_prepared_ = false;
    plus_group_warning_emitted_ = false;
  }
 
 private:
  Type PopulationSizeBinProb(const GrowthProducts<Type>& growth_products,
                             std::size_t year_index, std::size_t age_index,
                             std::size_t size_bin_index) const {
    if (!population_size_grid_ || population_size_grid_->n_bins == 0) {
      throw std::runtime_error(
          "GrowthDerivedSizeProvider requires a population biological size grid"
          " before computing ProbSize");
    }
 
    const Type mean_laa = growth_products.MeanLAA(year_index, age_index, 0);
    const Type sd_laa =
        fims_math::ad_max(growth_products.SdLAA(year_index, age_index, 0),
                          static_cast<Type>(1e-8));
 
    if (population_size_grid_->n_bins == 1) {
      return static_cast<Type>(1.0);
    }
 
    if (size_bin_index == 0) {
      const Type upper = static_cast<Type>(population_size_grid_->edges[1]);
      return fims_math::pnorm(upper, mean_laa, sd_laa);
    }
 
    if (size_bin_index + 1 == population_size_grid_->n_bins) {
      const Type lower = static_cast<Type>(
          population_size_grid_->edges[population_size_grid_->n_bins - 1]);
      return fims_math::pnorm(static_cast<Type>(2.0) * mean_laa - lower,
                              mean_laa, sd_laa);
    }
 
    const Type lower =
        static_cast<Type>(population_size_grid_->edges[size_bin_index]);
    const Type upper =
        static_cast<Type>(population_size_grid_->edges[size_bin_index + 1]);
 
    return fims_math::pnorm(upper, mean_laa, sd_laa) -
           fims_math::pnorm(lower, mean_laa, sd_laa);
  }
 
  const GrowthProducts<Type>& PreparedGrowthProducts() const {
    if (!growth_observation_) {
      throw std::runtime_error(
          "GrowthDerivedSizeProvider requires a linked growth object");
    }
 
    const GrowthProducts<Type>* growth_products =
        growth_observation_->TryGetPreparedGrowthProducts();
 
    if (!growth_products) {
      throw std::runtime_error(
          "GrowthDerivedSizeProvider requires prepared growth products");
    }
 
    return *growth_products;
  }
 
  std::shared_ptr<GrowthDerivedObservationBase<Type>> growth_observation_;
  const SizeGrid* population_size_grid_ = nullptr;
  std::size_t n_years_ = 0;
  std::size_t n_ages_ = 0;
  bool size_products_prepared_ = false;
  bool plus_group_warning_emitted_ = false;
  SizeProducts<Type> size_products_;
  fims::Vector<Type> mean_weight_at_age_;
};
 
}  // namespace fims_popdy
 
#endif /* FIMS_POPULATION_DYNAMICS_GROWTH_DERIVED_SIZE_PROVIDER_HPP */