From d021fc7bd600d6c529363611c69c7dfc6ef44cd5 Mon Sep 17 00:00:00 2001 From: Jim Ianelli Date: Fri, 24 Apr 2026 08:09:18 -0700 Subject: [PATCH] Add three-parameter double logistic selectivity --- NAMESPACE | 1 + R/FIMS-package.R | 1 + R/Rcpp_exports.R | 3 +- R/create_default_parameters.R | 25 ++- inst/WORDLIST | 3 + inst/include/common/fims_math.hpp | 40 ++++ .../include/interface/rcpp/rcpp_interface.hpp | 3 + .../rcpp/rcpp_objects/rcpp_selectivity.hpp | 211 ++++++++++++++++++ .../selectivity/functors/double_logistic3.hpp | 53 +++++ .../selectivity/selectivity.hpp | 1 + src/fims_modules.hpp | 12 + tests/gtest/CMakeLists.txt | 24 ++ ...c3_DoubleLogistic3Selectivity_Evaluate.cpp | 36 +++ .../gtest/test_fims_math_double_logistic3.cpp | 20 ++ tests/testthat/test-rcpp-selectivity.R | 42 ++++ 15 files changed, 472 insertions(+), 3 deletions(-) create mode 100644 inst/include/population_dynamics/selectivity/functors/double_logistic3.hpp create mode 100644 tests/gtest/test_DoubleLogistic3_DoubleLogistic3Selectivity_Evaluate.cpp create mode 100644 tests/gtest/test_fims_math_double_logistic3.cpp diff --git a/NAMESPACE b/NAMESPACE index f95145c36..513f381fd 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,6 +7,7 @@ export(CreateTMBModel) export(DlnormDistribution) export(DmultinomDistribution) export(DnormDistribution) +export(DoubleLogistic3Selectivity) export(DoubleLogisticSelectivity) export(EWAAGrowth) export(FIMSFit) diff --git a/R/FIMS-package.R b/R/FIMS-package.R index d52eec0d5..bc5ff4282 100644 --- a/R/FIMS-package.R +++ b/R/FIMS-package.R @@ -10,6 +10,7 @@ #' @export DmultinomDistribution #' @export DnormDistribution #' @export DoubleLogisticSelectivity +#' @export DoubleLogistic3Selectivity #' @export EWAAGrowth #' @export Fleet #' @export set_fixed diff --git a/R/Rcpp_exports.R b/R/Rcpp_exports.R index dfaf1d217..46a30205f 100644 --- a/R/Rcpp_exports.R +++ b/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) diff --git a/R/create_default_parameters.R b/R/create_default_parameters.R index b72cd2f67..0fbc0c7fb 100644 --- a/R/create_default_parameters.R +++ b/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" diff --git a/inst/WORDLIST b/inst/WORDLIST index c5a0f1cda..907268cdb 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -150,6 +150,9 @@ DoozyX DoubleAgeInput DoubleInput DoubleLogistic +DoubleLogistic3 +DoubleLogistic3Selectivity +DoubleLogistic3SelectivityInterface DoubleLogisticSel DoubleLogisticSelectivity DoubleLogisticSelectivityInterface diff --git a/inst/include/common/fims_math.hpp b/inst/include/common/fims_math.hpp index 9c38345ac..31c03bebf 100644 --- a/inst/include/common/fims_math.hpp +++ b/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(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 +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(2.0) * p1 + p2 + p3; + const Type log19 = fims_math::log(static_cast(19.0)); + const Type asc = + static_cast(1.0) / + (static_cast(1.0) + + exp(Type(-1.0) * log19 * (x - gamma1) / p1)); + const Type desc = + static_cast(1.0) - + static_cast(1.0) / + (static_cast(1.0) + + exp(Type(-1.0) * log19 * (x - gamma2) / p3)); + const Type normalized = asc * desc / + (static_cast(0.95) * static_cast(0.95)); + return fims_math::ad_min(normalized, static_cast(1.0), + static_cast(1e-12)); +} + /** * Returns the maximum between a and b in a continuous manner using: * diff --git a/inst/include/interface/rcpp/rcpp_interface.hpp b/inst/include/interface/rcpp/rcpp_interface.hpp index e35369c64..50fba5f6d 100644 --- a/inst/include/interface/rcpp/rcpp_interface.hpp +++ b/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(); diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp index 5656940ea..abc94cffe 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp +++ b/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(*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 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> info = + fims_info::Information::GetInstance(); + + fims_info::Information::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> sel = + std::dynamic_pointer_cast< + fims_popdy::DoubleLogistic3Selectivity>(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 + bool add_to_fims_tmb_internal() { + std::shared_ptr> info = + fims_info::Information::GetInstance(); + + std::shared_ptr> + selectivity = + std::make_shared>(); + + 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(); + this->add_to_fims_tmb_internal(); + + return true; + } + +#endif +}; + #endif diff --git a/inst/include/population_dynamics/selectivity/functors/double_logistic3.hpp b/inst/include/population_dynamics/selectivity/functors/double_logistic3.hpp new file mode 100644 index 000000000..7050579e9 --- /dev/null +++ b/inst/include/population_dynamics/selectivity/functors/double_logistic3.hpp @@ -0,0 +1,53 @@ +/** + * @file double_logistic3.hpp + * @brief Declares the DoubleLogistic3Selectivity class. + * @copyright This file is part of the NOAA, National Marine Fisheries Service + * Fisheries Integrated Modeling System project. See LICENSE in the source + * folder for reuse information. + */ +#ifndef POPULATION_DYNAMICS_SELECTIVITY_DOUBLE_LOGISTIC3_HPP +#define POPULATION_DYNAMICS_SELECTIVITY_DOUBLE_LOGISTIC3_HPP + +#include "../../../common/fims_math.hpp" +#include "../../../common/fims_vector.hpp" +#include "selectivity_base.hpp" + +namespace fims_popdy { + +/** + * @brief Three-parameter double logistic selectivity class. + */ +template +struct DoubleLogistic3Selectivity : public SelectivityBase { + fims::Vector p1; /**< Ascending limb width from 50% to 95%. */ + fims::Vector p2; /**< Horizontal shift of the ascending limb. */ + fims::Vector p3; /**< Descending limb width from 50% to 5%. */ + + DoubleLogistic3Selectivity() : SelectivityBase() {} + + virtual ~DoubleLogistic3Selectivity() {} + + /** + * @brief Evaluate the three-parameter double logistic selectivity function. + * + * @param x The independent variable, e.g., age or size. + */ + virtual const Type evaluate(const Type &x) { + return fims_math::double_logistic3(p1[0], p2[0], p3[0], x); + } + + /** + * @copydoc DoubleLogistic3Selectivity::evaluate(const Type &x) + * @param pos Position index, e.g., which year. If the index is out of bounds + * then it returns the first element, which is the time-invariant case. + */ + virtual const Type evaluate(const Type &x, size_t pos) { + return fims_math::double_logistic3( + p1.get_force_scalar(pos), p2.get_force_scalar(pos), + p3.get_force_scalar(pos), x); + } +}; + +} // namespace fims_popdy + +#endif /* POPULATION_DYNAMICS_SELECTIVITY_DOUBLE_LOGISTIC3_HPP */ diff --git a/inst/include/population_dynamics/selectivity/selectivity.hpp b/inst/include/population_dynamics/selectivity/selectivity.hpp index b5726d02b..d0c3568b5 100644 --- a/inst/include/population_dynamics/selectivity/selectivity.hpp +++ b/inst/include/population_dynamics/selectivity/selectivity.hpp @@ -12,6 +12,7 @@ #ifndef FIMS_POPULATION_DYNAMICS_SELECTIVITY_HPP #define FIMS_POPULATION_DYNAMICS_SELECTIVITY_HPP +#include "functors/double_logistic3.hpp" #include "functors/double_logistic.hpp" #include "functors/logistic.hpp" #include "functors/selectivity_base.hpp" diff --git a/src/fims_modules.hpp b/src/fims_modules.hpp index 77987dace..baf59b65a 100644 --- a/src/fims_modules.hpp +++ b/src/fims_modules.hpp @@ -396,6 +396,18 @@ RCPP_MODULE(fims) { .method("get_id", &DoubleLogisticSelectivityInterface::get_id) .method("evaluate", &DoubleLogisticSelectivityInterface::evaluate); + Rcpp::class_( + "DoubleLogistic3Selectivity", + "See " + "https://noaa-fims.github.io/FIMS/doxygen/" + "classDoubleLogistic3SelectivityInterface.html.") + .constructor() + .field("p1", &DoubleLogistic3SelectivityInterface::p1) + .field("p2", &DoubleLogistic3SelectivityInterface::p2) + .field("p3", &DoubleLogistic3SelectivityInterface::p3) + .method("get_id", &DoubleLogistic3SelectivityInterface::get_id) + .method("evaluate", &DoubleLogistic3SelectivityInterface::evaluate); + Rcpp::class_( "EWAAGrowth", "See " diff --git a/tests/gtest/CMakeLists.txt b/tests/gtest/CMakeLists.txt index 354e58bfe..e24dfe6e1 100644 --- a/tests/gtest/CMakeLists.txt +++ b/tests/gtest/CMakeLists.txt @@ -64,6 +64,18 @@ target_link_libraries(fims_math_double_logistic gtest_discover_tests(fims_math_double_logistic) +# test_fims_math_double_logistic3.cpp +add_executable(fims_math_double_logistic3 + test_fims_math_double_logistic3.cpp +) + +target_link_libraries(fims_math_double_logistic3 + gtest_main + fims_test +) + +gtest_discover_tests(fims_math_double_logistic3) + # test_fims_math_ad_fab_min_max.cpp add_executable(fims_math_fabs_min_max test_fims_math_ad_fabs_min_max.cpp @@ -290,6 +302,18 @@ target_link_libraries(DoubleLogistic_DoubleLogisticSelectivity_Evaluate gtest_discover_tests(DoubleLogistic_DoubleLogisticSelectivity_Evaluate) +# test_DoubleLogistic3_DoubleLogistic3Selectivity_Evaluate.cpp +add_executable(DoubleLogistic3_DoubleLogistic3Selectivity_Evaluate + test_DoubleLogistic3_DoubleLogistic3Selectivity_Evaluate.cpp +) + +target_link_libraries(DoubleLogistic3_DoubleLogistic3Selectivity_Evaluate + gtest_main + fims_test +) + +gtest_discover_tests(DoubleLogistic3_DoubleLogistic3Selectivity_Evaluate) + # test_Logistic_LogisticSelectivity_Evaluate.cpp add_executable(Logistic_LogisticSelectivity_Evaluate test_Logistic_LogisticSelectivity_Evaluate.cpp diff --git a/tests/gtest/test_DoubleLogistic3_DoubleLogistic3Selectivity_Evaluate.cpp b/tests/gtest/test_DoubleLogistic3_DoubleLogistic3Selectivity_Evaluate.cpp new file mode 100644 index 000000000..aa728f742 --- /dev/null +++ b/tests/gtest/test_DoubleLogistic3_DoubleLogistic3Selectivity_Evaluate.cpp @@ -0,0 +1,36 @@ +#include "gtest/gtest.h" +#include "population_dynamics/selectivity/functors/double_logistic3.hpp" + +namespace { + +TEST(DoubleLogistic3Selectivity_Evaluate, HandlesCorrectInput) { + fims_popdy::DoubleLogistic3Selectivity fishery_selectivity; + fishery_selectivity.p1.resize(1); + fishery_selectivity.p2.resize(1); + fishery_selectivity.p3.resize(1); + fishery_selectivity.p1[0] = 2.0; + fishery_selectivity.p2[0] = 4.0; + fishery_selectivity.p3[0] = 2.5; + + EXPECT_NEAR(fishery_selectivity.evaluate(9.0), 0.9350173, 0.0001); +} + +TEST(DoubleLogistic3Selectivity_Evaluate, HandlesThreeTimeSteps) { + fims_popdy::DoubleLogistic3Selectivity fishery_selectivity; + fishery_selectivity.p1.resize(3); + fishery_selectivity.p2.resize(1); + fishery_selectivity.p3.resize(1); + fishery_selectivity.p1[0] = 2.0; + fishery_selectivity.p1[1] = 1.0; + fishery_selectivity.p1[2] = 3.0; + fishery_selectivity.p2[0] = 4.0; + fishery_selectivity.p3[0] = 2.5; + + double expected_fishery[3] = {0.9350173, 0.3954430, 0.9560853}; + for (size_t pos = 0; pos < 3; ++pos) { + EXPECT_NEAR(fishery_selectivity.evaluate(9.0, pos), + expected_fishery[pos], 0.0001); + } +} + +} // namespace diff --git a/tests/gtest/test_fims_math_double_logistic3.cpp b/tests/gtest/test_fims_math_double_logistic3.cpp new file mode 100644 index 000000000..3e4dda67d --- /dev/null +++ b/tests/gtest/test_fims_math_double_logistic3.cpp @@ -0,0 +1,20 @@ +#include "common/fims_math.hpp" +#include "gtest/gtest.h" + +namespace { + +TEST(DoubleLogistic3, UseMultipleInputValues) { + double p1_value[3] = {2.0, 2.0, 1.0}; + double p2_value[3] = {4.0, 4.0, 4.0}; + double p3_value[3] = {2.5, 2.5, 2.5}; + double x_value[3] = {9.0, 5.0, 8.0}; + double expect_value[3] = {0.9350173, 0.2064477, 0.7124833}; + + for (size_t i = 0; i < 3; i++) { + EXPECT_NEAR(fims_math::double_logistic3(p1_value[i], p2_value[i], + p3_value[i], x_value[i]), + expect_value[i], 0.0001); + } +} + +} // namespace diff --git a/tests/testthat/test-rcpp-selectivity.R b/tests/testthat/test-rcpp-selectivity.R index ec3f449c7..8244a5a0c 100644 --- a/tests/testthat/test-rcpp-selectivity.R +++ b/tests/testthat/test-rcpp-selectivity.R @@ -99,6 +99,48 @@ test_that("rcpp double logistic selectivity works with correct inputs", { clear() }) +test_that("rcpp 3-parameter double logistic selectivity works with correct inputs", { + selectivity1 <- methods::new(DoubleLogistic3Selectivity) + + selectivity1$p1[1]$value <- 2.0 + selectivity1$p2[1]$value <- 4.0 + selectivity1$p3[1]$value <- 2.5 + selectivity1$p1[1]$estimation_type$set("fixed_effects") + + #' @description Test that `get_id()` for `DoubleLogistic3Selectivity` works. + expect_equal(selectivity1$get_id(), 1) + #' @description Test that the `p1` value is set to 2.0. + expect_equal(selectivity1$p1[1]$value, 2.0) + #' @description Test that the `p1` estimation type is set to "fixed_effects". + expect_equal(selectivity1$p1[1]$estimation_type$get(), "fixed_effects") + #' @description Test that `evaluate()` works for `DoubleLogistic3Selectivity`. + expect_equal( + selectivity1$evaluate(9.0), + 0.9350173, + tolerance = 0.0001 + ) + + selectivity2 <- methods::new(DoubleLogistic3Selectivity) + selectivity2$p1[1]$value <- 1.0 + selectivity2$p2[1]$value <- 4.0 + selectivity2$p3[1]$value <- 2.5 + selectivity2$p1[1]$estimation_type$set("random_effects") + selectivity2$p2[1]$estimation_type$set("random_effects") + selectivity2$p3[1]$estimation_type$set("random_effects") + + #' @description Test that `get_id()` for `DoubleLogistic3Selectivity` works when a second object is created. + expect_equal(selectivity2$get_id(), 2) + #' @description Test that the `p3` estimation type is set to "random_effects". + expect_equal(selectivity2$p3[1]$estimation_type$get(), "random_effects") + #' @description Test that `evaluate()` works for a second `DoubleLogistic3Selectivity`. + expect_equal( + selectivity2$evaluate(8.0), + 0.7124833, + tolerance = 0.0001 + ) + clear() +}) + ## Edge handling ---- test_that("rcpp selectivity returns correct outputs for edge cases", { # emptyLogistic