[codex] Add three-parameter double logistic selectivity

NAMESPACE

@@ -7,6 +7,7 @@ export(CreateTMBModel)
 export(DlnormDistribution)
 export(DmultinomDistribution)
 export(DnormDistribution)
+export(DoubleLogistic3Selectivity)
 export(DoubleLogisticSelectivity)
 export(EWAAGrowth)
 export(FIMSFit)

R/FIMS-package.R

@@ -10,6 +10,7 @@
 #' @export DmultinomDistribution
 #' @export DnormDistribution
 #' @export DoubleLogisticSelectivity
+#' @export DoubleLogistic3Selectivity
 #' @export EWAAGrowth
 #' @export Fleet
 #' @export set_fixed

R/Rcpp_exports.R

@@ -8,7 +8,7 @@
 #' [NOAA-FIMS C++ Documentation](https://noaa-fims.github.io/FIMS/doxygen/)
 #'
 #' @name Cpp_classes
-#' @aliases AgeComp BevertonHoltRecruitment CatchAtAge DlnormDistribution DmultinomDistribution DnormDistribution DoubleLogisticSelectivity EWAAGrowth Fleet Index Landings LengthComp LogDevsRecruitmentProcess LogRRecruitmentProcess LogisticMaturity LogisticSelectivity Parameter ParameterVector Population RealVector SharedInt SharedReal SharedString
+#' @aliases AgeComp BevertonHoltRecruitment CatchAtAge DlnormDistribution DmultinomDistribution DnormDistribution DoubleLogisticSelectivity DoubleLogistic3Selectivity EWAAGrowth Fleet Index Landings LengthComp LogDevsRecruitmentProcess LogRRecruitmentProcess LogisticMaturity LogisticSelectivity Parameter ParameterVector Population RealVector SharedInt SharedReal SharedString
 #'
 #' @details
 #' - [AgeComp](https://noaa-fims.github.io/FIMS/doxygen/classAgeCompDataInterface.html)
@@ -18,6 +18,7 @@
 #' - [DmultinomDistribution](https://noaa-fims.github.io/FIMS/doxygen/classDmultinomDistributionsInterface.html)
 #' - [DnormDistribution](https://noaa-fims.github.io/FIMS/doxygen/classDnormDistributionsInterface.html)
 #' - [DoubleLogisticSelectivity](https://noaa-fims.github.io/FIMS/doxygen/classDoubleLogisticSelectivityInterface.html)
+#' - [DoubleLogistic3Selectivity](https://noaa-fims.github.io/FIMS/doxygen/classDoubleLogistic3SelectivityInterface.html)
 #' - [EWAAGrowth](https://noaa-fims.github.io/FIMS/doxygen/classEWAAGrowthInterface.html)
 #' - [Fleet](https://noaa-fims.github.io/FIMS/doxygen/classFleetInterface.html)
 #' - [Index](https://noaa-fims.github.io/FIMS/doxygen/classIndexDataInterface.html)

R/create_default_parameters.R

@@ -375,6 +375,26 @@ create_default_DoubleLogistic <- function(module_name = NA_character_) {
     )
 }
 
+#' Create default 3-parameter double logistic parameters
+#'
+#' @description
+#' This function sets up default parameters for a 3-parameter double logistic
+#' function.
+#' @return
+#' A tibble containing default p1, p2, and p3 values and their estimation
+#' status.
+#' @noRd
+create_default_DoubleLogistic3 <- function(module_name = NA_character_) {
+  default <- create_default_parameters_template(n_parameters = 3) |>
+    dplyr::mutate(
+      module_name = !!module_name,
+      module_type = "DoubleLogistic3",
+      label = c("p1", "p2", "p3"),
+      value = c(2, 4, 2.5),
+      estimation_type = "fixed_effects"
+    )
+}
+
 #' Create default selectivity parameters
 #'
 #' @description
@@ -388,7 +408,7 @@ create_default_DoubleLogistic <- function(module_name = NA_character_) {
 #' of selectivity.
 #' @noRd
 create_default_selectivity <- function(
-  form = c("Logistic", "DoubleLogistic")
+  form = c("Logistic", "DoubleLogistic", "DoubleLogistic3")
 ) {
   # Input checks
   form <- rlang::arg_match(form)
@@ -397,7 +417,8 @@ create_default_selectivity <- function(
   # `switch`
   default <- switch(form,
     "Logistic" = create_default_Logistic(),
-    "DoubleLogistic" = create_default_DoubleLogistic()
+    "DoubleLogistic" = create_default_DoubleLogistic(),
+    "DoubleLogistic3" = create_default_DoubleLogistic3()
   ) |>
     dplyr::mutate(
       module_name = "Selectivity"

inst/WORDLIST

@@ -150,6 +150,9 @@ DoozyX
 DoubleAgeInput
 DoubleInput
 DoubleLogistic
+DoubleLogistic3
+DoubleLogistic3Selectivity
+DoubleLogistic3SelectivityInterface
 DoubleLogisticSel
 DoubleLogisticSelectivity
 DoubleLogisticSelectivityInterface

inst/include/common/fims_math.hpp

@@ -314,6 +314,46 @@ inline const Type ad_min(const Type &a, const Type &b, Type C = 1e-5) {
   return (a + b - fims_math::ad_fabs(a - b, C)) * static_cast<Type>(0.5);
 }
 
+/**
+ * @brief Three-parameter double logistic selectivity function.
+ *
+ * This form follows the parameterization:
+ * \f$ \gamma_1 = p_1 + p_2 \f$ and
+ * \f$ \gamma_2 = 2p_1 + p_2 + p_3 \f$.
+ *
+ * \f$ asc(x) = \frac{1}{1 + exp[-log(19)(x-\gamma_1)/p_1]} \f$
+ *
+ * \f$ desc(x) = 1 - \frac{1}{1 + exp[-log(19)(x-\gamma_2)/p_3]} \f$
+ *
+ * \f$ sel(x) = min(1, asc(x) desc(x) / 0.95^2) \f$
+ *
+ * @param p1 ascending limb width from 50% to 95% selectivity
+ * @param p2 horizontal shift of the ascending limb
+ * @param p3 descending limb width from 50% to 5% selectivity
+ * @param x the index the function should be evaluated at
+ * @return Selectivity at x, capped at one with a smooth minimum.
+ */
+template <class Type>
+inline const Type double_logistic3(const Type &p1, const Type &p2,
+                                   const Type &p3, const Type &x) {
+  const Type gamma1 = p1 + p2;
+  const Type gamma2 = static_cast<Type>(2.0) * p1 + p2 + p3;
+  const Type log19 = fims_math::log(static_cast<Type>(19.0));
+  const Type asc =
+      static_cast<Type>(1.0) /
+      (static_cast<Type>(1.0) +
+       exp(Type(-1.0) * log19 * (x - gamma1) / p1));
+  const Type desc =
+      static_cast<Type>(1.0) -
+      static_cast<Type>(1.0) /
+          (static_cast<Type>(1.0) +
+           exp(Type(-1.0) * log19 * (x - gamma2) / p3));
+  const Type normalized = asc * desc /
+                          (static_cast<Type>(0.95) * static_cast<Type>(0.95));
+  return fims_math::ad_min(normalized, static_cast<Type>(1.0),
+                           static_cast<Type>(1e-12));
+}
+
 /**
  * Returns the maximum between a and b in a continuous manner using:
  *

inst/include/interface/rcpp/rcpp_interface.hpp

@@ -309,6 +309,9 @@ void clear() {
   DoubleLogisticSelectivityInterface::id_g = 1;
   DoubleLogisticSelectivityInterface::live_objects.clear();
 
+  DoubleLogistic3SelectivityInterface::id_g = 1;
+  DoubleLogistic3SelectivityInterface::live_objects.clear();
+
   // rcpp_distribution.hpp
   DistributionsInterfaceBase::id_g = 1;
   DistributionsInterfaceBase::live_objects.clear();

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

@@ -606,4 +606,215 @@ class DoubleLogisticSelectivityInterface : public SelectivityInterfaceBase {
 #endif
 };
 
+/**
+ * @brief Rcpp interface for three-parameter double logistic selectivity.
+ */
+class DoubleLogistic3SelectivityInterface : public SelectivityInterfaceBase {
+ public:
+  ParameterVector p1; /**< Ascending limb width from 50% to 95%. */
+  ParameterVector p2; /**< Horizontal shift of the ascending limb. */
+  ParameterVector p3; /**< Descending limb width from 50% to 5%. */
+
+  DoubleLogistic3SelectivityInterface() : SelectivityInterfaceBase() {
+    SelectivityInterfaceBase::live_objects[this->id] =
+        std::make_shared<DoubleLogistic3SelectivityInterface>(*this);
+    FIMSRcppInterfaceBase::fims_interface_objects.push_back(
+        SelectivityInterfaceBase::live_objects[this->id]);
+  }
+
+  DoubleLogistic3SelectivityInterface(
+      const DoubleLogistic3SelectivityInterface &other)
+      : SelectivityInterfaceBase(other),
+        p1(other.p1),
+        p2(other.p2),
+        p3(other.p3) {}
+
+  virtual ~DoubleLogistic3SelectivityInterface() {}
+
+  virtual uint32_t get_id() { return this->id; }
+
+  virtual double evaluate(double x) {
+    fims_popdy::DoubleLogistic3Selectivity<double> DoubleLogistic3Sel;
+    DoubleLogistic3Sel.p1.resize(1);
+    DoubleLogistic3Sel.p1[0] = this->p1[0].initial_value_m;
+    DoubleLogistic3Sel.p2.resize(1);
+    DoubleLogistic3Sel.p2[0] = this->p2[0].initial_value_m;
+    DoubleLogistic3Sel.p3.resize(1);
+    DoubleLogistic3Sel.p3[0] = this->p3[0].initial_value_m;
+    return DoubleLogistic3Sel.evaluate(x);
+  }
+
+  virtual void finalize() {
+    if (this->finalized) {
+      FIMS_WARNING_LOG("Double Logistic 3 Selectivity  " +
+                       fims::to_string(this->id) +
+                       " has been finalized already.");
+    }
+
+    this->finalized = true;
+
+    std::shared_ptr<fims_info::Information<double>> info =
+        fims_info::Information<double>::GetInstance();
+
+    fims_info::Information<double>::selectivity_models_iterator it =
+        info->selectivity_models.find(this->id);
+    if (it == info->selectivity_models.end()) {
+      FIMS_WARNING_LOG("Double Logistic 3 Selectivity " +
+                       fims::to_string(this->id) +
+                       " not found in Information.");
+      return;
+    } else {
+      std::shared_ptr<fims_popdy::DoubleLogistic3Selectivity<double>> sel =
+          std::dynamic_pointer_cast<
+              fims_popdy::DoubleLogistic3Selectivity<double>>(it->second);
+
+      for (size_t i = 0; i < p1.size(); i++) {
+        if (this->p1[i].estimation_type_m.get() == "constant") {
+          this->p1[i].final_value_m = this->p1[i].initial_value_m;
+        } else {
+          this->p1[i].final_value_m = sel->p1[i];
+        }
+      }
+
+      for (size_t i = 0; i < p2.size(); i++) {
+        if (this->p2[i].estimation_type_m.get() == "constant") {
+          this->p2[i].final_value_m = this->p2[i].initial_value_m;
+        } else {
+          this->p2[i].final_value_m = sel->p2[i];
+        }
+      }
+
+      for (size_t i = 0; i < p3.size(); i++) {
+        if (this->p3[i].estimation_type_m.get() == "constant") {
+          this->p3[i].final_value_m = this->p3[i].initial_value_m;
+        } else {
+          this->p3[i].final_value_m = sel->p3[i];
+        }
+      }
+    }
+  }
+
+  virtual std::string to_json() {
+    std::stringstream ss;
+
+    ss << "{\n";
+    ss << " \"module_name\": \"Selectivity\",\n";
+    ss << " \"module_type\": \"DoubleLogistic3\",\n";
+    ss << " \"module_id\": " << this->id << ",\n";
+
+    ss << " \"parameters\":[\n{\n";
+    ss << "   \"name\": \"p1\",\n";
+    ss << "   \"id\":" << this->p1.id_m << ",\n";
+    ss << "   \"type\": \"vector\",\n";
+    ss << " \"dimensionality\": {\n";
+    ss << "  \"header\": [null],\n";
+    ss << "  \"dimensions\": [1]\n},\n";
+    ss << "   \"values\":" << this->p1 << "},\n";
+
+    ss << "{\n";
+    ss << "   \"name\": \"p2\",\n";
+    ss << "   \"id\":" << this->p2.id_m << ",\n";
+    ss << "   \"type\": \"vector\",\n";
+    ss << " \"dimensionality\": {\n";
+    ss << "  \"header\": [null],\n";
+    ss << "  \"dimensions\": [1]\n},\n";
+    ss << "   \"values\":" << this->p2 << "},\n";
+
+    ss << "{\n";
+    ss << "   \"name\": \"p3\",\n";
+    ss << "   \"id\":" << this->p3.id_m << ",\n";
+    ss << "   \"type\": \"vector\",\n";
+    ss << " \"dimensionality\": {\n";
+    ss << "  \"header\": [null],\n";
+    ss << "  \"dimensions\": [1]\n},\n";
+    ss << "   \"values\":" << this->p3 << "}]\n";
+
+    ss << "}";
+
+    return ss.str();
+  }
+
+#ifdef TMB_MODEL
+
+  template <typename Type>
+  bool add_to_fims_tmb_internal() {
+    std::shared_ptr<fims_info::Information<Type>> info =
+        fims_info::Information<Type>::GetInstance();
+
+    std::shared_ptr<fims_popdy::DoubleLogistic3Selectivity<Type>>
+        selectivity =
+            std::make_shared<fims_popdy::DoubleLogistic3Selectivity<Type>>();
+
+    std::stringstream ss;
+    selectivity->id = this->id;
+
+    selectivity->p1.resize(this->p1.size());
+    for (size_t i = 0; i < this->p1.size(); i++) {
+      selectivity->p1[i] = this->p1[i].initial_value_m;
+      if (this->p1[i].estimation_type_m.get() == "fixed_effects") {
+        ss.str("");
+        ss << "Selectivity." << this->id << ".p1." << this->p1[i].id_m;
+        info->RegisterParameterName(ss.str());
+        info->RegisterParameter(selectivity->p1[i]);
+      }
+      if (this->p1[i].estimation_type_m.get() == "random_effects") {
+        ss.str("");
+        ss << "Selectivity." << this->id << ".p1." << this->p1[i].id_m;
+        info->RegisterRandomEffectName(ss.str());
+        info->RegisterRandomEffect(selectivity->p1[i]);
+      }
+    }
+    info->variable_map[this->p1.id_m] = &(selectivity)->p1;
+
+    selectivity->p2.resize(this->p2.size());
+    for (size_t i = 0; i < this->p2.size(); i++) {
+      selectivity->p2[i] = this->p2[i].initial_value_m;
+      if (this->p2[i].estimation_type_m.get() == "fixed_effects") {
+        ss.str("");
+        ss << "Selectivity." << this->id << ".p2." << this->p2[i].id_m;
+        info->RegisterParameterName(ss.str());
+        info->RegisterParameter(selectivity->p2[i]);
+      }
+      if (this->p2[i].estimation_type_m.get() == "random_effects") {
+        ss.str("");
+        ss << "Selectivity." << this->id << ".p2." << this->p2[i].id_m;
+        info->RegisterRandomEffectName(ss.str());
+        info->RegisterRandomEffect(selectivity->p2[i]);
+      }
+    }
+    info->variable_map[this->p2.id_m] = &(selectivity)->p2;
+
+    selectivity->p3.resize(this->p3.size());
+    for (size_t i = 0; i < this->p3.size(); i++) {
+      selectivity->p3[i] = this->p3[i].initial_value_m;
+      if (this->p3[i].estimation_type_m.get() == "fixed_effects") {
+        ss.str("");
+        ss << "Selectivity." << this->id << ".p3." << this->p3[i].id_m;
+        info->RegisterParameterName(ss.str());
+        info->RegisterParameter(selectivity->p3[i]);
+      }
+      if (this->p3[i].estimation_type_m.get() == "random_effects") {
+        ss.str("");
+        ss << "Selectivity." << this->id << ".p3." << this->p3[i].id_m;
+        info->RegisterRandomEffectName(ss.str());
+        info->RegisterRandomEffect(selectivity->p3[i]);
+      }
+    }
+    info->variable_map[this->p3.id_m] = &(selectivity)->p3;
+
+    info->selectivity_models[selectivity->id] = selectivity;
+
+    return true;
+  }
+
+  virtual bool add_to_fims_tmb() {
+    this->add_to_fims_tmb_internal<TMB_FIMS_REAL_TYPE>();
+    this->add_to_fims_tmb_internal<TMBAD_FIMS_TYPE>();
+
+    return true;
+  }
+
+#endif
+};
+
 #endif