growth_model_adapter.hpp

Go to the documentation of this file.

#ifndef POPULATION_DYNAMICS_GROWTH_MODEL_ADAPTER_HPP
#define POPULATION_DYNAMICS_GROWTH_MODEL_ADAPTER_HPP
 
#include <cmath>
#include <memory>
#include <stdexcept>
 
#include "../../common/fims_vector.hpp"
#include "../../common/fims_math.hpp"
#include "growth_model.hpp"
#include "functors/growth_base.hpp"
#include "functors/vonb_schnute.hpp"
 
namespace fims_popdy {

template <typename Type>
class GrowthDerivedObservationBase : public GrowthBase<Type> {
 public:
  GrowthDerivedObservationBase() : GrowthBase<Type>() {}
  virtual ~GrowthDerivedObservationBase() = default;

  virtual void SetAgeOffset(double min_age) = 0;

  virtual void Initialize(std::size_t n_years, std::size_t n_ages,
                          std::size_t n_sexes = 1) = 0;

  virtual bool SupportsAgeToLengthConversionDerived() const = 0;

  virtual void PrepareGrowthProducts() = 0;

  virtual const GrowthProducts<Type>* TryGetPreparedGrowthProducts() const = 0;

  virtual Type EvaluateWeightAtLength(const Type& length) const = 0;
};

template <typename Type>
class VonBertalanffySchnuteGrowthModelAdapter
    : public GrowthDerivedObservationBase<Type> {
 public:
  VonBertalanffySchnuteGrowthModelAdapter()
      : GrowthDerivedObservationBase<Type>() {}

  enum class LengthReferenceParameterization {
    kEstimatedReferenceLengths = 0,
    kConstantLengthYoungEstimatedLengthOld,
    kEstimatedLengthYoungBelowConstantLengthOld,
    kBothConstant
  };

  void UseEstimatedReferenceMeanLengths() {
    length_reference_parameterization_ =
        LengthReferenceParameterization::kEstimatedReferenceLengths;
    growth_products_prepared_ = false;
  }

  void UseConstantMeanLengthYoungWithEstimatedMeanLengthOld(
      Type constant_mean_length_young) {
    length_reference_parameterization_ =
        LengthReferenceParameterization::kConstantLengthYoungEstimatedLengthOld;
    constant_mean_length_young_ = constant_mean_length_young;
    growth_products_prepared_ = false;
  }

  void UseEstimatedMeanLengthYoungBelowConstantMeanLengthOld(
      Type constant_mean_length_old) {
    length_reference_parameterization_ = LengthReferenceParameterization::
        kEstimatedLengthYoungBelowConstantLengthOld;
    constant_mean_length_old_ = constant_mean_length_old;
    growth_products_prepared_ = false;
  }

  void UseConstantReferenceMeanLengths(Type constant_mean_length_young,
                                       Type constant_mean_length_old) {
    length_reference_parameterization_ =
        LengthReferenceParameterization::kBothConstant;
    constant_mean_length_young_ = constant_mean_length_young;
    constant_mean_length_old_ = constant_mean_length_old;
    growth_products_prepared_ = false;
  }

  fims::Vector<Type>& MeanLengthYoungVector() {
    use_param_vectors_ = true;
    vb_params_set_ = true;
    growth_products_prepared_ = false;
    return mean_length_young_vector_;
  }

  fims::Vector<Type>& MeanLengthOldVector() {
    use_param_vectors_ = true;
    vb_params_set_ = true;
    growth_products_prepared_ = false;
    return mean_length_old_vector_;
  }

  fims::Vector<Type>& GrowthCoefficientVector() {
    use_param_vectors_ = true;
    vb_params_set_ = true;
    growth_products_prepared_ = false;
    return growth_coefficient_vector_;
  }

  fims::Vector<Type>& ReferenceAgeForLengthYoungVector() {
    use_param_vectors_ = true;
    vb_params_set_ = true;
    growth_products_prepared_ = false;
    return reference_age_for_length_young_vector_;
  }

  fims::Vector<Type>& ReferenceAgeForLengthOldVector() {
    use_param_vectors_ = true;
    vb_params_set_ = true;
    growth_products_prepared_ = false;
    return reference_age_for_length_old_vector_;
  }

  fims::Vector<Type>& LengthWeightAVector() {
    use_param_vectors_ = true;
    lw_params_set_ = true;
    growth_products_prepared_ = false;
    return length_weight_a_vector_;
  }

  fims::Vector<Type>& LengthWeightBVector() {
    use_param_vectors_ = true;
    lw_params_set_ = true;
    growth_products_prepared_ = false;
    return length_weight_b_vector_;
  }

  fims::Vector<Type>& LengthAtAgeSdYoungVector() {
    use_param_vectors_ = true;
    growth_products_prepared_ = false;
    return length_at_age_sd_young_vector_;
  }

  fims::Vector<Type>& LengthAtAgeSdOldVector() {
    use_param_vectors_ = true;
    growth_products_prepared_ = false;
    return length_at_age_sd_old_vector_;
  }

  fims::Vector<Type>& LogSdLengthAtRefAgeYoungVector() {
    use_param_vectors_ = true;
    growth_products_prepared_ = false;
    return log_sd_mean_length_young_vector_;
  }

  fims::Vector<Type>& LogSdLengthAtRefAgeOldVector() {
    use_param_vectors_ = true;
    growth_products_prepared_ = false;
    return log_sd_mean_length_old_vector_;
  }

  fims::Vector<Type>& LogSdGrowthCoefficientVector() {
    use_param_vectors_ = true;
    growth_products_prepared_ = false;
    return log_sd_growth_coefficient_vector_;
  }

  fims::Vector<Type>& LogitCorrLengthAtRefAgeYoungLengthAtRefAgeOldVector() {
    use_param_vectors_ = true;
    growth_products_prepared_ = false;
    return mean_length_young_mean_length_old_logit_corr_vector_;
  }

  fims::Vector<Type>& LogitCorrLengthAtRefAgeYoungKVector() {
    use_param_vectors_ = true;
    growth_products_prepared_ = false;
    return mean_length_young_growth_coefficient_logit_corr_vector_;
  }

  fims::Vector<Type>& LogitCorrLengthAtRefAgeOldKVector() {
    use_param_vectors_ = true;
    growth_products_prepared_ = false;
    return mean_length_old_growth_coefficient_logit_corr_vector_;
  }

  void SetAgeOffset(double min_age) override {
    age_offset_ = min_age;
    age_offset_set_ = true;
    growth_products_prepared_ = false;
    if (model_) {
      model_->SetAgeOffset(static_cast<Type>(age_offset_));
    }
  }

  void Initialize(std::size_t n_years, std::size_t n_ages,
                  std::size_t n_sexes = 1) override {
    EnsureParamsSet();
    if (n_sexes != 1) {
      throw std::runtime_error(
          "VonBertalanffySchnuteGrowthModelAdapter currently supports n_sexes "
          "== 1");
    }
    n_years_ = n_years;
    n_ages_ = n_ages;
    n_sexes_ = n_sexes;
    growth_products_prepared_ = false;
    model_ = std::make_shared<GrowthModel<Type>>(n_years, n_ages, n_sexes);
    SyncParamsToModel();
    if (age_offset_set_) {
      model_->SetAgeOffset(static_cast<Type>(age_offset_));
    }
  }

  bool SupportsAgeToLengthConversionDerived() const override { return true; }
 
  virtual const Type evaluate(int year, const double& a) override {
    if (a < 0.0) {
      throw std::runtime_error("Negative age not supported");
    }
    const double a_round = std::round(a);
    const double tol = 1e-8;
    if (std::fabs(a - a_round) > tol) {
      throw std::runtime_error("Non-integer age not supported yet");
    }
    EnsureParamsSet();
    const Type ref_age_young = CurrentReferenceAgeForLengthYoung();
    const Type ref_age_old = CurrentReferenceAgeForLengthOld();
    if (ref_age_old <= ref_age_young) {
      throw std::runtime_error(
          "VonBertalanffySchnuteGrowth reference_age_for_length_old must be > "
          "reference_age_for_length_young");
    }
 
    if (!model_) {
      return EvaluateWithFunctor(year, a);
    }
 
    SyncParamsToModel();
    model_->Prepare();
    const auto& p = model_->GetProducts();
    const double offset = age_offset_set_ ? age_offset_ : 0.0;
    const double age_index_raw = a_round - offset;
    const double age_index_round = std::round(age_index_raw);
    if (std::fabs(age_index_raw - age_index_round) <= tol &&
        age_index_round >= 0.0) {
      const std::size_t age_index = static_cast<std::size_t>(age_index_round);
      if (age_index < p.n_ages) {
        const std::size_t year_index =
            (year >= 0 && static_cast<std::size_t>(year) < p.n_years)
                ? static_cast<std::size_t>(year)
                : static_cast<std::size_t>(0);
        return p.MeanWAA(year_index, age_index, 0);
      }
    }
 
    // Outside cached model age bins: fall back to direct functor evaluation.
    return EvaluateWithFunctor(year, a);
  }

  void PrepareGrowthProducts() override {
    growth_products_prepared_ = false;
    if (!model_) {
      if (n_ages_ == 0) {
        throw std::runtime_error("Growth model not initialized; n_ages is 0");
      }
      Initialize(n_years_ == 0 ? 1 : n_years_, n_ages_,
                 n_sexes_ == 0 ? 1 : n_sexes_);
    }
    SyncParamsToModel();
    model_->Prepare();
    growth_products_prepared_ = true;
  }

  const GrowthProducts<Type>* TryGetPreparedGrowthProducts() const override {
    if (!model_ || !growth_products_prepared_) {
      return nullptr;
    }
    return &(model_->GetProducts());
  }

  Type EvaluateWeightAtLength(const Type& length) const override {
    EnsureParamsSet();
    const Type length_safe = fims_math::ad_max(length, static_cast<Type>(1e-8));
    return CurrentLengthWeightA() *
           fims_math::pow(length_safe, CurrentLengthWeightB());
  }
 
 private:
  // Stored parameter vectors on their working scales. Positive parameters
  // use log scale, reference ages stay on the natural scale, and
  // correlation terms use transformed working-scale values.
  fims::Vector<Type> mean_length_young_vector_;
  fims::Vector<Type> mean_length_old_vector_;
  fims::Vector<Type> growth_coefficient_vector_;
  fims::Vector<Type> reference_age_for_length_young_vector_;
  fims::Vector<Type> reference_age_for_length_old_vector_;
  fims::Vector<Type> length_weight_a_vector_;
  fims::Vector<Type> length_weight_b_vector_;
  fims::Vector<Type> length_at_age_sd_young_vector_;
  fims::Vector<Type> length_at_age_sd_old_vector_;
  fims::Vector<Type> log_sd_mean_length_young_vector_;
  fims::Vector<Type> log_sd_mean_length_old_vector_;
  fims::Vector<Type> log_sd_growth_coefficient_vector_;
  fims::Vector<Type> mean_length_young_mean_length_old_logit_corr_vector_;
  fims::Vector<Type> mean_length_young_growth_coefficient_logit_corr_vector_;
  fims::Vector<Type> mean_length_old_growth_coefficient_logit_corr_vector_;
  bool use_param_vectors_ = false;
  LengthReferenceParameterization length_reference_parameterization_ =
      LengthReferenceParameterization::kEstimatedReferenceLengths;
  Type constant_mean_length_young_ = Type(0.0);
  Type constant_mean_length_old_ = Type(0.0);
  std::size_t n_years_ = 0;
  std::size_t n_ages_ = 0;
  std::size_t n_sexes_ = 1;
  double age_offset_ = 0.0;
  bool age_offset_set_ = false;
  bool vb_params_set_ = false;
  bool lw_params_set_ = false;
  bool growth_products_prepared_ = false;
 
  mutable std::shared_ptr<GrowthModel<Type>> model_;
 
  Type EvaluateWithFunctor(int year, const double& a) const {
    EnsureParamsSet();
    fims_popdy::VonBertalanffySchnuteGrowth<Type> vb;
    vb.mean_length_young = CurrentMeanLengthYoung();
    vb.mean_length_old = CurrentMeanLengthOld();
    vb.growth_coefficient = CurrentGrowthCoefficient();
    vb.reference_age_for_length_young = CurrentReferenceAgeForLengthYoung();
    vb.reference_age_for_length_old = CurrentReferenceAgeForLengthOld();
    vb.length_weight_a = CurrentLengthWeightA();
    vb.length_weight_b = CurrentLengthWeightB();
    return vb.evaluate(year, a);
  }
 
  void SyncParamsToModel() const {
    if (!model_) return;
    EnsureParamsSet();
 
    model_->SetVonBertalanffySchnuteParameters(
        CurrentMeanLengthYoung(), CurrentMeanLengthOld(),
        CurrentGrowthCoefficient(), CurrentReferenceAgeForLengthYoung(),
        CurrentReferenceAgeForLengthOld());
    model_->SetLengthWeightParameters(CurrentLengthWeightA(),
                                      CurrentLengthWeightB());
 
    if (HasInterpolationSdInputs()) {
      model_->SetLengthSdParams(CurrentLengthAtAgeSdAtReferenceAge1(),
                                CurrentLengthAtAgeSdAtReferenceAge2());
    }
 
    if (HasStructuredDeltaMethodInputs()) {
      model_->SetGrowthParameterCovariance(
          CurrentMeanLengthYoungVariance(),
          CurrentMeanLengthYoungLengthAtRefAgeOldCovariance(),
          CurrentMeanLengthYoungKCovariance(), CurrentMeanLengthOldVariance(),
          CurrentMeanLengthOldKCovariance(),
          CurrentGrowthCoefficientVariance());
    } else {
      model_->ClearGrowthParameterCovariance();
    }
  }
 
  Type CurrentMeanLengthYoungStorage() const {
    return fims_math::exp(mean_length_young_vector_[0]);
  }
 
  Type CurrentMeanLengthOldStorage() const {
    return fims_math::exp(mean_length_old_vector_[0]);
  }
 
  Type CurrentMeanLengthYoung() const {
    switch (length_reference_parameterization_) {
      case LengthReferenceParameterization::kEstimatedReferenceLengths:
        return CurrentMeanLengthYoungStorage();
 
      case LengthReferenceParameterization::
          kConstantLengthYoungEstimatedLengthOld:
        return constant_mean_length_young_;
 
      case LengthReferenceParameterization::
          kEstimatedLengthYoungBelowConstantLengthOld:
        return fims_math::inv_logit(static_cast<Type>(0.0),
                                    constant_mean_length_old_,
                                    mean_length_young_vector_[0]);
 
      case LengthReferenceParameterization::kBothConstant:
        return constant_mean_length_young_;
    }
 
    return CurrentMeanLengthYoungStorage();
  }
 
  Type CurrentMeanLengthOld() const {
    switch (length_reference_parameterization_) {
      case LengthReferenceParameterization::kEstimatedReferenceLengths:
      case LengthReferenceParameterization::
          kConstantLengthYoungEstimatedLengthOld:
        return CurrentMeanLengthOldStorage();
 
      case LengthReferenceParameterization::
          kEstimatedLengthYoungBelowConstantLengthOld:
        return constant_mean_length_old_;
 
      case LengthReferenceParameterization::kBothConstant:
        return constant_mean_length_old_;
    }
 
    return CurrentMeanLengthOldStorage();
  }
 
  Type CurrentGrowthCoefficient() const {
    return fims_math::exp(growth_coefficient_vector_[0]);
  }
  Type CurrentReferenceAgeForLengthYoung() const {
    return reference_age_for_length_young_vector_[0];
  }
  Type CurrentReferenceAgeForLengthOld() const {
    return reference_age_for_length_old_vector_[0];
  }
  Type CurrentLengthWeightA() const {
    return fims_math::exp(length_weight_a_vector_[0]);
  }
  Type CurrentLengthWeightB() const {
    return fims_math::exp(length_weight_b_vector_[0]);
  }
  Type CurrentLengthAtAgeSdAtReferenceAge1() const {
    return fims_math::exp(length_at_age_sd_young_vector_[0]);
  }
  Type CurrentLengthAtAgeSdAtReferenceAge2() const {
    return fims_math::exp(length_at_age_sd_old_vector_[0]);
  }
  Type CurrentSdLengthAtRefAgeYoung() const {
    return fims_math::exp(log_sd_mean_length_young_vector_[0]);
  }
 
  Type CurrentSdLengthAtRefAgeOld() const {
    return fims_math::exp(log_sd_mean_length_old_vector_[0]);
  }
 
  Type CurrentSdGrowthCoefficient() const {
    return fims_math::exp(log_sd_growth_coefficient_vector_[0]);
  }
 
  Type CurrentCorrLengthAtRefAgeYoungLengthAtRefAgeOld() const {
    return fims_math::inv_logit(
        static_cast<Type>(-1.0), static_cast<Type>(1.0),
        mean_length_young_mean_length_old_logit_corr_vector_[0]);
  }
 
  Type CurrentCorrLengthAtRefAgeYoungK() const {
    return fims_math::inv_logit(
        static_cast<Type>(-1.0), static_cast<Type>(1.0),
        mean_length_young_growth_coefficient_logit_corr_vector_[0]);
  }
 
  Type CurrentCorrLengthAtRefAgeOldK() const {
    return fims_math::inv_logit(
        static_cast<Type>(-1.0), static_cast<Type>(1.0),
        mean_length_old_growth_coefficient_logit_corr_vector_[0]);
  }
 
  Type CurrentMeanLengthYoungVariance() const {
    const Type sd = CurrentSdLengthAtRefAgeYoung();
    return sd * sd;
  }
 
  Type CurrentMeanLengthOldVariance() const {
    const Type sd = CurrentSdLengthAtRefAgeOld();
    return sd * sd;
  }
 
  Type CurrentGrowthCoefficientVariance() const {
    const Type sd = CurrentSdGrowthCoefficient();
    return sd * sd;
  }
 
  Type CurrentMeanLengthYoungLengthAtRefAgeOldCovariance() const {
    return CurrentCorrLengthAtRefAgeYoungLengthAtRefAgeOld() *
           CurrentSdLengthAtRefAgeYoung() * CurrentSdLengthAtRefAgeOld();
  }
 
  Type CurrentMeanLengthYoungKCovariance() const {
    return CurrentCorrLengthAtRefAgeYoungK() * CurrentSdLengthAtRefAgeYoung() *
           CurrentSdGrowthCoefficient();
  }
 
  Type CurrentMeanLengthOldKCovariance() const {
    return CurrentCorrLengthAtRefAgeOldK() * CurrentSdLengthAtRefAgeOld() *
           CurrentSdGrowthCoefficient();
  }
 
  bool HasInterpolationSdInputs() const {
    return length_at_age_sd_young_vector_.size() > 0 ||
           length_at_age_sd_old_vector_.size() > 0;
  }
 
  bool HasAnyStructuredDeltaMethodInput() const {
    return log_sd_mean_length_young_vector_.size() > 0 ||
           log_sd_mean_length_old_vector_.size() > 0 ||
           log_sd_growth_coefficient_vector_.size() > 0 ||
           mean_length_young_mean_length_old_logit_corr_vector_.size() > 0 ||
           mean_length_young_growth_coefficient_logit_corr_vector_.size() > 0 ||
           mean_length_old_growth_coefficient_logit_corr_vector_.size() > 0;
  }
 
  bool HasStructuredDeltaMethodInputs() const {
    return log_sd_mean_length_young_vector_.size() > 0 &&
           log_sd_mean_length_old_vector_.size() > 0 &&
           log_sd_growth_coefficient_vector_.size() > 0 &&
           mean_length_young_mean_length_old_logit_corr_vector_.size() > 0 &&
           mean_length_young_growth_coefficient_logit_corr_vector_.size() > 0 &&
           mean_length_old_growth_coefficient_logit_corr_vector_.size() > 0;
  }
 
  void EnsureParamsSet() const {
    if (!use_param_vectors_ || mean_length_young_vector_.size() < 1 ||
        mean_length_old_vector_.size() < 1 ||
        growth_coefficient_vector_.size() < 1 ||
        reference_age_for_length_young_vector_.size() < 1 ||
        reference_age_for_length_old_vector_.size() < 1 ||
        length_weight_a_vector_.size() < 1 ||
        length_weight_b_vector_.size() < 1) {
      throw std::runtime_error(
          "VonBertalanffySchnuteGrowth parameters not set");
    }
 
    if (mean_length_young_vector_.size() != 1 ||
        mean_length_old_vector_.size() != 1 ||
        growth_coefficient_vector_.size() != 1 ||
        reference_age_for_length_young_vector_.size() != 1 ||
        reference_age_for_length_old_vector_.size() != 1 ||
        length_weight_a_vector_.size() != 1 ||
        length_weight_b_vector_.size() != 1) {
      throw std::runtime_error(
          "VonBertalanffySchnuteGrowthModelAdapter currently supports a single "
          "growth pattern; expected size 1 for von Bertalanffy--Schnute and "
          "length-weight parameter vectors");
    }
 
    const bool has_sd = HasInterpolationSdInputs();
    const bool has_any_structured_delta = HasAnyStructuredDeltaMethodInput();
    const bool has_structured_delta = HasStructuredDeltaMethodInputs();
 
    if (!has_sd && !has_any_structured_delta) {
      throw std::runtime_error(
          "VonBertalanffySchnuteGrowthModelAdapter requires either "
          "both length_at_age_sd_young and length_at_age_sd_old, "
          "or the structured delta-method growth variability inputs");
    }
 
    if (has_any_structured_delta && !has_structured_delta) {
      throw std::runtime_error(
          "VonBertalanffySchnuteGrowthModelAdapter requires all six structured "
          "delta-method variability inputs when using that path");
    }
 
    if (has_sd && has_structured_delta) {
      throw std::runtime_error(
          "von Bertalanffy--Schnute growth adapter requires variability inputs "
          "for exactly one supported path. Supply either the interpolation "
          "inputs length_at_age_sd_young and length_at_age_sd_old or the full "
          "delta-method variability inputs, but not both");
    }
 
    if (has_sd && (length_at_age_sd_young_vector_.size() != 1 ||
                   length_at_age_sd_old_vector_.size() != 1)) {
      throw std::runtime_error(
          "VonBertalanffySchnuteGrowthModelAdapter requires exactly one "
          "length_at_age_sd_young value and one length_at_age_sd_old value");
    }
 
    if (has_structured_delta &&
        (log_sd_mean_length_young_vector_.size() != 1 ||
         log_sd_mean_length_old_vector_.size() != 1 ||
         log_sd_growth_coefficient_vector_.size() != 1 ||
         mean_length_young_mean_length_old_logit_corr_vector_.size() != 1 ||
         mean_length_young_growth_coefficient_logit_corr_vector_.size() != 1 ||
         mean_length_old_growth_coefficient_logit_corr_vector_.size() != 1)) {
      throw std::runtime_error(
          "VonBertalanffySchnuteGrowthModelAdapter currently supports a single "
          "structured delta-method variability parameter set; expected "
          "size 1 for each structured uncertainty input");
    }
  }
};
 
}  // namespace fims_popdy
 
#endif  // POPULATION_DYNAMICS_GROWTH_MODEL_ADAPTER_HPP