From b88bae0884a9d81344089df882a5fd3c5c3a2dcb Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA <85530309+Andrea-Havron-NOAA@users.noreply.github.com> Date: Sun, 12 Apr 2026 16:55:56 +0000 Subject: [PATCH 1/4] Add generic transformation functionality C++ changes - create transformation enum struct and distribution enum struct in common/types.hpp - create common/fims_transformations.hpp with transformation functions - ApplyTransformation, ApplyBackTransformation, TransformPrior, AddLogJacobian - restructure variable_map to hold transformation metadata for parameters - add AddLogJacobian call to model.hpp and create do_mcmc flag - add Prepare function to distributions_base to calculate prior transform for each iteration Interface Changes - add new functions: add_prior to ParameterVector and add_shared_prior to rcpp_interface.hpp - add code to parse R formula and pass transformation information and prior informationinto the backend - add set_variable_map function to pass values and transformation metadata from Rcpp to C++ - modify calls to variable_map and set input transformations if parameters are not logged Tests - add C++ tests for fims_transformation.hpp - update R tests and vignettes cleanup code - remove min/max field and functions --- CMakeLists.txt | 9 +- NAMESPACE | 1 + R/FIMS-package.R | 1 + R/Rcpp_exports.R | 1 + inst/include/common/fims_transformations.hpp | 340 ++++++++++++++++++ inst/include/common/fims_vector.hpp | 1 + inst/include/common/information.hpp | 60 +++- inst/include/common/model.hpp | 15 + inst/include/common/types.hpp | 105 ++++++ .../functors/density_components_base.hpp | 43 ++- .../distributions/functors/normal_lpdf.hpp | 2 +- .../include/interface/rcpp/rcpp_interface.hpp | 229 +++++++++++- .../rcpp/rcpp_objects/rcpp_distribution.hpp | 15 +- .../rcpp/rcpp_objects/rcpp_fleet.hpp | 4 +- .../rcpp/rcpp_objects/rcpp_interface_base.hpp | 64 +++- .../rcpp/rcpp_objects/rcpp_models.hpp | 12 +- .../rcpp/rcpp_objects/rcpp_population.hpp | 8 +- .../rcpp/rcpp_objects/rcpp_recruitment.hpp | 17 +- .../rcpp/rcpp_objects/rcpp_selectivity.hpp | 12 +- man/Cpp_functions.Rd | 1 + src/FIMS.cpp | 3 + src/fims_modules.hpp | 4 + src/init.hpp | 1 + tests/gtest/CMakeLists.txt | 24 ++ ...ms_transformations_ApplyTransformation.cpp | 244 +++++++++++++ ...st_fims_transformations_TransformPrior.cpp | 134 +++++++ tests/gtest/test_info_setup_priors.cpp | 14 +- .../gtest/test_info_setup_random_effects.cpp | 4 +- .../helper-integration-tests-setup-function.R | 1 + ...tegration-fims-bayesian-prior-predictive.R | 16 +- tests/testthat/test-rcpp-distribution.R | 1 + 31 files changed, 1317 insertions(+), 69 deletions(-) create mode 100644 inst/include/common/fims_transformations.hpp create mode 100644 inst/include/common/types.hpp create mode 100644 tests/gtest/test_fims_transformations_ApplyTransformation.cpp create mode 100644 tests/gtest/test_fims_transformations_TransformPrior.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 36f8a81ed..2a8a17079 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -121,12 +121,19 @@ enable_testing() # Creates an interface library fims_test add_library(fims_test INTERFACE) +# Get the exact include path directly from R +EXECUTE_PROCESS( + COMMAND R "--slave" "--no-save" "-e" "cat(R.home('include'))" + OUTPUT_VARIABLE R_INCLUDE_DIR + OUTPUT_STRIP_TRAILING_WHITESPACE +) + # Add include directories to the fims_test target target_include_directories(fims_test INTERFACE inst/include ${rcpp_SOURCE_DIR}/inst/include - ${R_HOME}/include + ${R_INCLUDE_DIR} ) # Add compile definition STD_LIB to the fims_test target. diff --git a/NAMESPACE b/NAMESPACE index f95145c36..663d0ffea 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -26,6 +26,7 @@ export(RealVector) export(SharedInt) export(SharedReal) export(SharedString) +export(add_shared_prior) export(clear) export(create_default_configurations) export(create_default_parameters) diff --git a/R/FIMS-package.R b/R/FIMS-package.R index d52eec0d5..507e8cc9f 100644 --- a/R/FIMS-package.R +++ b/R/FIMS-package.R @@ -12,6 +12,7 @@ #' @export DoubleLogisticSelectivity #' @export EWAAGrowth #' @export Fleet +#' @export add_shared_prior #' @export set_fixed #' @export get_fixed #' @export set_random diff --git a/R/Rcpp_exports.R b/R/Rcpp_exports.R index dfaf1d217..be98ae892 100644 --- a/R/Rcpp_exports.R +++ b/R/Rcpp_exports.R @@ -50,6 +50,7 @@ NULL #' #' @details #' - [clear](https://noaa-fims.github.io/FIMS/doxygen/rcpp__interface_8hpp.html) +#' - [add_shared_prior](https://noaa-fims.github.io/FIMS/doxygen/add__shared__prior_8hpp.html) #' - [get_fixed](https://noaa-fims.github.io/FIMS/doxygen/rcpp__interface_8hpp.html) #' - [get_log](https://noaa-fims.github.io/FIMS/doxygen/rcpp__interface_8hpp.html) #' - [get_log_errors](https://noaa-fims.github.io/FIMS/doxygen/rcpp__interface_8hpp.html) diff --git a/inst/include/common/fims_transformations.hpp b/inst/include/common/fims_transformations.hpp new file mode 100644 index 000000000..8fe276f7a --- /dev/null +++ b/inst/include/common/fims_transformations.hpp @@ -0,0 +1,340 @@ +/** + * @file fims_transformations.hpp + * @brief Defines transformations for parameters in FIMS. + + * + * @details This file provides functions for applying forward and inverse + * transformations to parameters, transforming parameters between input + * and prior spaces, and computing log Jacobian adjustments for + * change-of-variables corrections in MCMC sampling. + * + * The transformation system supports the following transformations: + * - `identity`: No transformation, parameter is on the natural scale. + * - `exp`: Exponential transformation, parameter is on the log scale. + * - `log`: Log transformation, parameter is on the natural scale. + * - `logit`: Logit transformation with optional bounds. + * - `square`: Square transformation. + * - `sqrt`: Square root transformation. + * + * @see fims::Transformation + * @see fims_math.hpp + */ + +#ifndef FIMS_COMMON_FIMS_TRANSFORMATIONS_HPP +#define FIMS_COMMON_FIMS_TRANSFORMATIONS_HPP + +#include "fims_math.hpp" +#include "types.hpp" + +namespace fims_transformations { + + /** + * @brief Applies a forward transformation to a scalar value. + * + * @details Maps a value from natural scale to the transformed scale + * specified by the transformation argument. The forward transformations + * and are: + * + * | Label | Forward | + * |----------|------------------| + * | identity | x | + * | log | log(x) | + * | exp | exp(x) | + * | logit | logit(x, lo, hi) | + * | square | x^2 | + * | sqrt | sqrt(x) | + * + * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). + * @param value The input value on the natural scale. + * @param transformation A fims::Transformation struct specifying the + * transformation label and any required arguments (e.g. logit bounds). + * @return The transformed value. + * @throws std::invalid_argument if the transformation label is not supported. + */ + template + inline Type ApplyTransformation( + const Type& value, + const fims::Transformation& transformation) + { + + const auto label = transformation.label; + const auto args = transformation.args; + + switch(label) { + case fims::Transformation::Label::identity: + return value; + case fims::Transformation::Label::exp: + return exp(value); + case fims::Transformation::Label::log: + return log(value); + case fims::Transformation::Label::logit: + return fims_math::logit( + Type(args.lower), + Type(args.upper), value); + case fims::Transformation::Label::square: + return value * value; + case fims::Transformation::Label::sqrt: + return sqrt(value); + // Add more cases as needed + default: + throw std::invalid_argument("Unknown transformation label"); + } + } + + + /** + * @brief Applies an inverse (back) transformation to a scalar value. + * + * @details Maps a value from the transformed scale back to the natural + * scale. This is the inverse of ApplyTransformation(). For example, if + * the input transformation is `log`, then ApplyBackTransformation applies + * `exp` to recover the natural-scale value. + * + * | Label | Back Transformation | + * |----------|----------------------| + * | identity | x | + * | log | exp(x) | + * | exp | log(x) | + * | logit | inv_logit(x, lo, hi) | + * | square | sqrt(x) | + * | sqrt | x^2 | + * + * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). + * @param value The input value on the transformed scale. + * @param transformation A fims::Transformation struct specifying the + * transformation label and any required arguments (e.g. logit bounds). + * @return The back-transformed value on the natural scale. + * @throws std::invalid_argument if the transformation label is not supported. + */ + template + inline Type ApplyBackTransformation( + const Type& value, + const fims::Transformation& transformation) { + + const auto label = transformation.label; + const auto args = transformation.args; + + Type transformed_value; + switch (label) { + case fims::Transformation::Label::identity: + transformed_value = value; + break; + case fims::Transformation::Label::exp: + transformed_value = fims_math::log(value); + break; + case fims::Transformation::Label::log: + transformed_value = fims_math::exp(value); + break; + case fims::Transformation::Label::logit: + transformed_value = fims_math::inv_logit(Type(args.lower), + Type(args.upper), value); + break; + case fims::Transformation::Label::square: + transformed_value = fims_math::sqrt(value); + break; + case fims::Transformation::Label::sqrt: + transformed_value = value * value; + break; + default: + throw std::invalid_argument( + std::string("Unknown transformation applied to a parameter, ") + + TransformationLabelToString(label) + + std::string(". Valid transformations are identity, exp, log, logit, square, and sqrt.")); + break; + } + + return transformed_value; + } + + + /** + * @brief Transforms a vector of parameters from input space to prior space. + * + * @details Applies the composite transformation: + * input space -> natural scale -> prior space + * + * This is used when evaluating a prior distribution where the parameter + * is estimated in one space (e.g. log scale) but the prior is defined in + * a different space (e.g. natural scale or variance scale). + * + * If input and prior transformations are identical, the input vector is + * returned unchanged. + * + * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). + * @param input_value A vector of parameter values in input space. + * @param input The transformation applied to the parameter in input space + * (e.g. log if the parameter is stored as log(sd)). + * @param prior The transformation applied to the parameter in prior space + * (e.g. square if the prior is on variance = sd^2). + * @return A vector of parameter values transformed to prior space. + */ + template + inline fims::Vector TransformPrior(const fims::Vector& input_value, + fims::Transformation input, + fims::Transformation prior) { + + if (input.label == prior.label) { + return input_value; + } else { + + size_t n = input_value.size(); + + fims::Vector natural_parameter; + fims::Vector prior_parameter; + natural_parameter.resize(n); + prior_parameter.resize(n); + + for(size_t i = 0; i < n; i ++) { + natural_parameter[i] = + ApplyBackTransformation(input_value[i], input); + prior_parameter[i] = + ApplyTransformation(natural_parameter[i], prior); + } + + return prior_parameter; + } + } + + + /** + * @brief Transforms a scalar parameter from input space to prior space. + * + * @details Scalar overload of TransformPrior() for use in get_observed() + * where a single element is needed. Applies the composite transformation: + * input space -> natural scale -> prior space + * + * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). + * @param input_value A scalar parameter value in input space. + * @param input The transformation applied to the parameter in input space. + * @param prior The transformation applied to the parameter in prior space. + * @return The scalar parameter value transformed to prior space. + */ + template + inline Type TransformPrior(const Type& input_value, + const fims::Transformation& input, + const fims::Transformation& prior) { + + Type natural_parameter = ApplyBackTransformation(input_value, input); + return ApplyTransformation(natural_parameter, + prior); + } + + /** + * @brief Computes the log absolute determinant of the Jacobian for a + * change-of-variables correction in MCMC sampling. + * + * @details When a parameter is estimated in one space (input space) but + * the prior is placed in a different space (prior space), MCMC sampling + * requires a change-of-variables correction. This function computes + * log |det J| where J is the Jacobian of the composite transformation + * from input space to prior space. + * + * The Jacobian is computed using TMBad's sparse Jacobian functionality + * via SpJacFun(), which avoids allocating a dense n x n matrix. For + * element-wise transformations (the common case), the Jacobian is + * diagonal and only n non-zero entries are stored. For future multivariate + * or simplex transformations, the sparse pattern is detected automatically + * during AD taping. + * + * @note This function is only available when compiled with TMB + * (TMB_MODEL defined). For non-TMB builds it throws std::invalid_argument. + * + * @note Stan handles the Jacobian for the transformation from unconstrained + * sampling space to bounded input space automatically. This function only + * handles the additional correction from input space to prior space. + * + * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). + * @param input_value A vector of parameter values in input space. + * @param input The transformation applied to the parameter in input space + * (e.g. log if the parameter is stored as log(sd)). + * @param prior The transformation applied to the parameter in prior space + * (e.g. square if the prior is on variance = sd^2). + * @return The log absolute determinant of the Jacobian, log |det J|. + * @throws std::invalid_argument if called outside a TMB model context. + * + * @example + * // Parameter estimated as log(sd), prior on variance (sd^2) + * // Jacobian adjustment: log(2) + 2*log(sd) + * fims::Transformation log_trans, square_trans; + * log_trans.label = fims::Transformation::Label::log; + * square_trans.label = fims::Transformation::Label::square; + * Type jac = AddLogJacobian(log_sd_values, log_trans, square_trans); + */ + template + inline Type AddLogJacobian(const fims::Vector& input_value, + fims::Transformation input, + fims::Transformation prior) { + + #ifdef TMB_MODEL + + size_t n = input_value.size(); + // Step 1: Extract double values for taping + std::vector x_val(n); + for (size_t i = 0; i < n; i++) x_val[i] = asDouble(input_value[i]); + + // Step 2: Tape the composite transformation into an ADFun + TMBad::ADFun<> f; + { + // Start recording + + std::vector x_ad(n); + for (size_t i = 0; i < n; i++) x_ad[i] = x_val[i]; + + f.glob.ad_start(); + f.glob.Independent(x_ad); // vector version, not scalar + + std::vector y_ad(n); + for (size_t i = 0; i < n; i++) { + TMBad::ad_aug natural = ApplyBackTransformation(x_ad[i], input); + y_ad[i] = ApplyTransformation(natural, prior); + } + + TMBad::Dependent(y_ad); // vector version + f.glob.ad_stop(); + } + + // Step 3: Generate sparse Jacobian function + TMBad::Sparse> spjac = f.SpJacFun(); + + // Step 4: Evaluate sparse Jacobian at current values + std::vector jac_vals = spjac(x_val); + + // Step 5: Compute log abs det from sparse values + // For diagonal case, all non-zeros are diagonal entries + // For general case, need to use spjac.i and spjac.j for row/col indices + Type log_abs_det = Type(0.0); + + if (n == spjac.i.size()) { + // Diagonal case - all non-zeros are on diagonal + for (size_t k = 0; k < jac_vals.size(); k++) { + log_abs_det += fims_math::log(fims_math::ad_fabs(Type(jac_vals[k]))); + } + } else { + // General case - need log determinant of sparse matrix + // Build Eigen sparse matrix from spjac pattern + Eigen::SparseMatrix J(n, n); + std::vector> triplets; + for (size_t k = 0; k < jac_vals.size(); k++) { + triplets.push_back(Eigen::Triplet( + spjac.i[k], spjac.j[k], jac_vals[k])); + } + J.setFromTriplets(triplets.begin(), triplets.end()); + + // For sparse log determinant use Eigen's sparse LU + Eigen::SparseLU> solver; + solver.compute(J); + log_abs_det = Type(solver.logAbsDeterminant()); + } + + return log_abs_det; + + #else + + throw std::invalid_argument( + std::string("Jacobian adjustments currently only work for TMB models.") + ); + #endif //TMB_MODEL + } + +} // namespace fims_transformations +#endif /* FIMS_COMMON_FIMS_TRANSFORMATIONS_HPP */ \ No newline at end of file diff --git a/inst/include/common/fims_vector.hpp b/inst/include/common/fims_vector.hpp index 176ae7c11..6b899afd3 100644 --- a/inst/include/common/fims_vector.hpp +++ b/inst/include/common/fims_vector.hpp @@ -8,6 +8,7 @@ #ifndef FIMS_VECTOR_HPP #define FIMS_VECTOR_HPP +#include "types.hpp" #include "../interface/interface.hpp" #include #include diff --git a/inst/include/common/information.hpp b/inst/include/common/information.hpp index 55bfb2c29..fd6f43107 100644 --- a/inst/include/common/information.hpp +++ b/inst/include/common/information.hpp @@ -138,11 +138,55 @@ class Information { typedef typename std::unordered_map< uint32_t, std::shared_ptr>>::iterator model_map_iterator; /**< iterator for variable map>*/ + +/** + * @brief A structure to hold a pointer to a parameter vector and its transformation + * metadata for use in the variable map. + * + * @details Each entry in the variable map corresponds to a single parameter + * vector and stores a pointer to the parameter values along with two + * transformation labels: + * - `input_transformation`: the transformation applied to the parameter + * in the input space (e.g. log, logit). This is the space in which + * the parameter is estimated. + * - `prior_transformation`: the transformation applied to the parameter + * in the prior space (e.g. identity, square). This is the space in + * which the prior distribution is defined. + */ +struct VariableMapEntry { + fims::Vector* variable = nullptr; + fims::Transformation input_transformation; + fims::Transformation prior_transformation; + +/** + * @brief Constructor for VariableMapEntry. + */ +VariableMapEntry() { + input_transformation.label = fims::Transformation::Label::log; + prior_transformation.label = fims::Transformation::Label::log; +} + +/** + * @brief Constructor for VariableMapEntry with all fields initialized. + * + * @param variable Pointer to the fims::Vector holding the parameter values. + * @param input_transformation The transformation applied to the parameter + * in the input space (e.g. log, logit). + * @param prior_transformation The transformation applied to the parameter + * in the prior space (e.g. identity, square). + */ + VariableMapEntry(fims::Vector* variable, + fims::Transformation input_transformation, + fims::Transformation prior_transformation) + : variable(variable), + input_transformation(input_transformation), + prior_transformation(prior_transformation) {} +}; - std::unordered_map*> + std::unordered_map variable_map; /***>::iterator + typedef typename std::unordered_map::iterator variable_map_iterator; /**< iterator for variable map>*/ Information() {} @@ -294,12 +338,16 @@ class Information { FIMS_INFO_LOG("Link prior from distribution " + fims::to_string(d->id) + " to parameter " + fims::to_string(d->key[0])); d->priors.resize(d->key.size()); + d->input_transformation.resize(d->key.size()); + d->prior_transformation.resize(d->key.size()); for (size_t i = 0; i < d->key.size(); i++) { FIMS_INFO_LOG("Link prior from distribution " + fims::to_string(d->id) + " to parameter " + fims::to_string(d->key[0])); vmit = this->variable_map.find(d->key[i]); - d->priors[i] = (*vmit).second; + d->priors[i] = (*vmit).second.variable; + d->input_transformation[i] = &(*vmit).second.input_transformation; + d->prior_transformation[i] = &(*vmit).second.prior_transformation; } FIMS_INFO_LOG("Prior size for distribution " + fims::to_string(d->id) + "is: " + fims::to_string(d->observed_values.size())); @@ -324,10 +372,10 @@ class Information { fims::to_string(d->id) + " to derived value " + fims::to_string(d->key[0])); vmit = this->variable_map.find(d->key[0]); - d->re = (*vmit).second; + d->re = (*vmit).second.variable; if (d->key.size() == 2) { vmit = this->variable_map.find(d->key[1]); - d->re_expected_values = (*vmit).second; + d->re_expected_values = (*vmit).second.variable; } else { d->re_expected_values = &d->expected_values; } @@ -355,7 +403,7 @@ class Information { fims::to_string(d->id) + " to derived value " + fims::to_string(d->key[0])); vmit = this->variable_map.find(d->key[0]); - d->data_expected_values = (*vmit).second; + d->data_expected_values = (*vmit).second.variable; FIMS_INFO_LOG( "Expected value size for distribution " + fims::to_string(d->id) + " is: " + fims::to_string((*d->data_expected_values).size())); diff --git a/inst/include/common/model.hpp b/inst/include/common/model.hpp index 998752691..43d5ea1de 100644 --- a/inst/include/common/model.hpp +++ b/inst/include/common/model.hpp @@ -13,6 +13,7 @@ #include #include "information.hpp" +#include "fims_transformations.hpp" namespace fims_model { @@ -27,6 +28,7 @@ class Model { // may need singleton std::shared_ptr> fims_information; /**< Create a shared fims_information as a pointer to Information*/ + bool jacobian_flag; /** * @brief Construct a new Model object. @@ -99,8 +101,21 @@ class Model { // may need singleton #ifdef TMB_MODEL d->of = this->of; #endif + d->Prepare(); if (d->input_type == "prior") { nll_vec[nll_vec_idx] = -d->evaluate(); + if (this->jacobian_flag) { + for(size_t i=0; i < d->key.size(); i++) { + if ((*(d->input_transformation[i])).label != + (*(d->prior_transformation[i])).label) { + nll_vec[nll_vec_idx] += + -fims_transformations::AddLogJacobian( + *(d->priors[i]), + *(d->input_transformation[i]), + *(d->prior_transformation[i])); + } + } + } jnll += nll_vec[nll_vec_idx]; n_priors += 1; nll_vec_idx += 1; diff --git a/inst/include/common/types.hpp b/inst/include/common/types.hpp new file mode 100644 index 000000000..737851db5 --- /dev/null +++ b/inst/include/common/types.hpp @@ -0,0 +1,105 @@ +/** + * @file types.hpp + * @brief Centralized definition of core model types, enumerations, and structures. + * @details This file serves as a lightweight, dependency-free foundation for the + * fims namespace. + * @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 FIMS_COMMON_TYPES_HPP +#define FIMS_COMMON_TYPES_HPP + +namespace fims { + +/** + * @brief A data structure for transforming parameters + * + * @details The `Transformation` struct defines the type of transformation to + * apply to a parameter and any necessary arguments for that transformation. + * This allows for flexible parameter transformations within the FIMS framework. + */ +struct Transformation { + /** + * @brief An enumeration of supported transformation types. + * @details Supported transformations include: + * - `identity`: No transformation, parameter is on the natural scale. + * - `exp`: Exponential transformation, parameter is on the multiplicative scale. + * - `log`: Log transformation, parameter is on the log scale. + * - `logit`: Logit transformation, parameter is on the logit scale. + * - `square`: Square transformation, parameter is on the squared scale. + * - `sqrt`: Square root transformation, parameter is on the square root scale. + */ + enum class Label { + identity, + exp, + log, + logit, + square, + sqrt + }; + /** + * @brief A structure to hold arguments for transformations that require them. + * @details For the logit transformation, `lower` and `upper` specify the + * bounds of the parameter in natural space. + */ + struct Args { + //used by logit + double lower = 0.0; + double upper = 1.0; + }; + + + Label label = Label::log; + Args args{}; + + }; + + static const char* TransformationLabelToString( + fims::Transformation::Label label) { + switch (label) { + case fims::Transformation::Label::identity: return "identity"; + case fims::Transformation::Label::exp: return "exp"; + case fims::Transformation::Label::log: return "log"; + case fims::Transformation::Label::logit: return "logit"; + case fims::Transformation::Label::square: return "square"; + case fims::Transformation::Label::sqrt: return "sqrt"; + default: return "unknown"; + } + } + + /** + * @brief A structure to hold distribution type information for priors. + * @details Supported distributions include: + * - `Normal`: Normal distribution. + * - `Lognormal`: Lognormal distribution. + * - `Gamma`: Gamma distribution. + * - `InvGamma`: Inverse gamma distribution. + * - `Multinom`: Multinomial distribution. + */ +struct Distribution { + /** + * @brief An enumeration of supported distribution types. + */ + enum class Label { + Normal, + Lognormal, + Gamma, + InvGamma + }; + + Label label = Label::Normal; +}; + +fims::Distribution::Label StringToDistributionLabel(const std::string& name) { + if (name == "Normal") return fims::Distribution::Label::Normal; + if (name == "Lognormal") return fims::Distribution::Label::Lognormal; + if (name == "Gamma") return fims::Distribution::Label::Gamma; + if (name == "InvGamma") return fims::Distribution::Label::InvGamma; + throw std::invalid_argument("Unsupported distribution: " + name); +} + +} // namespace fims + +#endif /* FIMS_COMMON_TYPES_HPP */ \ No newline at end of file diff --git a/inst/include/distributions/functors/density_components_base.hpp b/inst/include/distributions/functors/density_components_base.hpp index 1e052a1d6..5bda02799 100644 --- a/inst/include/distributions/functors/density_components_base.hpp +++ b/inst/include/distributions/functors/density_components_base.hpp @@ -18,8 +18,10 @@ #include "../../common/data_object.hpp" #include "../../common/model_object.hpp" #include "../../interface/interface.hpp" +#include "../../common/fims_transformations.hpp" #include "../../common/fims_vector.hpp" #include "../../common/fims_math.hpp" +#include "../../common/types.hpp" namespace fims_distributions { @@ -55,6 +57,15 @@ struct DensityComponentBase : public fims_model_object::FIMSObject { /** @brief Vector of pointers where each entry points to a prior parameter. */ std::vector*> priors; + /** @brief Pointer to input transformation */ + std::vector input_transformation; + + /** @brief Pointer to prior transformation */ + std::vector prior_transformation; + + /** @brief Prior values after transformation */ + fims::Vector transformed_priors; + /** * @brief Input value of distribution function for priors or random effects. */ @@ -71,8 +82,28 @@ struct DensityComponentBase : public fims_model_object::FIMSObject { */ std::string use_mean = fims::to_string("no"); - // std::shared_ptr> expected; - // // Expected value of distribution function. + virtual void Prepare() { + if (this->input_type == "prior") { + if (priors.size() == 0) { + throw std::runtime_error("No priors defined for this distribution."); + } else if (priors.size() == 1) { + this->transformed_priors.resize((*priors[0]).size()); + this->transformed_priors = + fims_transformations::TransformPrior(*(priors[0]), + *(input_transformation[0]), + *(prior_transformation[0])); + } else if (priors.size() > 1) { + size_t n = priors.size(); + this->transformed_priors.resize(n); + for (size_t i = 0; i < n; i++) { + transformed_priors[i] = + fims_transformations::TransformPrior((*(priors[i]))[0], + *(input_transformation[i]), + *(prior_transformation[i])); + } + } + } +} /** * @brief Retrieve one observed value based on `input_type`. @@ -88,13 +119,7 @@ struct DensityComponentBase : public fims_model_object::FIMSObject { return (*re)[i]; } if (this->input_type == "prior") { - if (priors.size() == 0) { - throw std::runtime_error("No priors defined for this distribution."); - } else if (priors.size() == 1) { - return (*(priors[0]))[i]; - } else if (priors.size() > 1) { - return (*(priors[i]))[0]; - } + return transformed_priors[i]; } return observed_values[i]; } diff --git a/inst/include/distributions/functors/normal_lpdf.hpp b/inst/include/distributions/functors/normal_lpdf.hpp index a0b60aa91..24238ddfa 100644 --- a/inst/include/distributions/functors/normal_lpdf.hpp +++ b/inst/include/distributions/functors/normal_lpdf.hpp @@ -114,7 +114,7 @@ struct NormalLPDF : public DensityComponentBase { } if (this->input_type == "random_effects") { (*this->re)[i] = rnorm(this->get_expected(i), - fims_math::exp(log_sd.get_force_scalar(i))); + fims_math::exp(log_sd.get_force_scalar(i))); } if (this->input_type == "prior") { (*(this->priors[i]))[0] = diff --git a/inst/include/interface/rcpp/rcpp_interface.hpp b/inst/include/interface/rcpp/rcpp_interface.hpp index e35369c64..993fbae69 100644 --- a/inst/include/interface/rcpp/rcpp_interface.hpp +++ b/inst/include/interface/rcpp/rcpp_interface.hpp @@ -8,6 +8,8 @@ #ifndef FIMS_INTERFACE_RCPP_INTERFACE_HPP #define FIMS_INTERFACE_RCPP_INTERFACE_HPP #include "../../common/model.hpp" +#include "../../common/fims_math.hpp" +#include "../../common/types.hpp" #include "../../utilities/fims_json.hpp" #include "rcpp_objects/rcpp_data.hpp" #include "rcpp_objects/rcpp_distribution.hpp" @@ -94,11 +96,11 @@ bool CreateTMBModel() { info0->CheckModel(); info->CreateModel(); - +/* // instantiate the model? TODO: Ask Matthew what this does std::shared_ptr> m0 = fims_model::Model::GetInstance(); - +*/ return true; } @@ -226,6 +228,229 @@ Rcpp::List get_random_names(Rcpp::List pars) { return pars; } +/** + * @brief Classifies an R formula's Left-Hand Side (LHS) and extracts the transformation. + * * @details This helper function inspects the Abstract Syntax Tree (AST) of the raw + * R expression pointer (SEXP) passed from the formula's left-hand side. It maps + * recognized R mathematical operations (like `log()`, `exp()`, or exponentiation `^2`) + * to their respective domain-specific strongly typed enums defined in the fims namespace. + * * Supported AST structural mappings include: + * - \b Symbol (e.g., \code{sd}): Maps to \c Label::identity + * - \b Language Call (e.g., \code{log(sd)}): Maps to \c Label::log + * - \b Language Call (e.g., \code{exp(sd)}): Maps to \c Label::exp + * - \b Language Call (e.g., \code{logit(sd)}): Maps to \c Label::logit + * - \b Language Call (e.g., \code{sqrt(sd)}): Maps to \c Label::sqrt + * - \b Infix Exponent (e.g., \code{sd^2}): Maps to \c Label::square + * + * @param lhs_raw A raw \code{SEXP} representing the left-hand side node of an R formula. + * * @throws std::runtime_error (via \code{Rcpp::stop}) if the expression tree structure + * is unrecognizable or contains an unmapped transformation function. + * * @return A \code{fims::Transformation::Label} with the enumeration token. + */ +fims::Transformation classify_and_extract_transformation(SEXP lhs_raw) { + + fims::Transformation trans_config; + + if (Rcpp::is(lhs_raw)) { + trans_config.label = fims::Transformation::Label::identity; + return trans_config; + } + + if (Rcpp::is(lhs_raw)) { + Rcpp::Language lhs_lang = Rcpp::as(lhs_raw); + Rcpp::CharacterVector op_name = Rcpp::as(Rcpp::as(lhs_lang[0])); + + // Match R operators/functions to your fims::Transformation::Label enums + if (op_name[0] == "log") { + trans_config.label = fims::Transformation::Label::log; + return trans_config; + } + else if (op_name[0] == "exp") { + trans_config.label = fims::Transformation::Label::exp; + return trans_config; + } + else if (op_name[0] == "logit") { + trans_config.label = fims::Transformation::Label::logit; + // Extract bounds if provided, otherwise default to 0 and 1 + if (lhs_lang.size() >= 3) { + trans_config.args.lower = Rcpp::as(lhs_lang[2]); + } + if (lhs_lang.size() >= 4) { + trans_config.args.upper = Rcpp::as(lhs_lang[3]); + } + // Validate bounds + if (trans_config.args.lower >= trans_config.args.upper) { + Rcpp::stop("Logit transformation requires lower < upper."); + } + return trans_config; + } + else if (op_name[0] == "sqrt") { + trans_config.label = fims::Transformation::Label::sqrt; + return trans_config; + } + else if (op_name[0] == "^") { + int power = Rcpp::as(lhs_lang[2]); + if (power == 2) { + trans_config.label = fims::Transformation::Label::square; + return trans_config; + } + } + // Add reciprocal parsing here if you still want to route 1/x variants + } + + // Fallback/Error state if formula layout isn't supported + Rcpp::stop("Unsupported Left-Hand Side (LHS) transformation operator."); +} + + +/** + * @brief A structure to hold the parsed components of a distribution formula. + */ +struct FormulaComponents { + fims::Transformation transformation; + fims::Distribution::Label distribution; + std::vector hyperparameters; +}; + +//' Parse a Distributional Formula +//' +//' Parses an R formula specifying a target variable and its prior/likelihood distribution +//' along with its parameters (e.g., \code{y ~ dnorm(0, 1)}). +//' +//' @param f A standard R \code{Formula} object. It must follow the structure +//' \code{variable ~ distribution(param1, param2, ...)}. +//' +//' @details +//' The function unpacks the R formula by treating it as an abstract syntax tree (AST) +//' via the \code{Rcpp::Language} class: +//' \itemize{ +//' \item \strong{Operator (\code{[0]}):} The tilde (\code{~}) operator. +//' \item \strong{LHS (\code{[1]}):} Extracted as a \code{Symbol} and converted to a character vector representing the target variable name. +//' \item \strong{RHS (\code{[2]}):} Treated as a nested \code{Language} call where the head (\code{[0]}) is the distribution name, and subsequent elements are the numeric parameters. +//' } +//' +//' @return A named \code{Rcpp::List} containing three elements: +//' \itemize{ +//' \item \code{variable}: A character vector holding the name of the LHS variable. +//' \item \code{distribution}: A character vector holding the name of the RHS distribution function. +//' \item \code{hyperparameters}: A numeric vector containing the extracted hyperparameter values. +//' } +FormulaComponents parse_distribution_formula(Rcpp::Formula f) { + + // Convert the Formula to a standard standard R language object (Call) + Rcpp::Language formula = Rcpp::as(f); + + // Parse the formula + // R formulas are structured as a tree where element 0 is the operator `~` + // element 1 is the Left-Hand Side (LHS), and element 2 is the Right-Hand Side (RHS) + + SEXP lhs_raw = formula[1]; + + // Extract the transformation configuration and variable string + fims::Transformation transform = classify_and_extract_transformation(lhs_raw); + + // 2. Parse Right-Hand Side (RHS) + if (!Rcpp::is(formula[2])) { + Rcpp::stop("Right-Hand Side (RHS) must be a distribution function call."); + } + Rcpp::Language rhs = Rcpp::as(formula[2]); + + // The first element of a function call language object is the function name itself + // Convert formula->Symbol->CharacterVector->enum + Rcpp::Symbol distribution_symbol = Rcpp::as(rhs[0]); + Rcpp::CharacterVector distribution_string = Rcpp::as(distribution_symbol); + + // 3. Extract the numeric parameters from the function arguments + int num_args = rhs.size() - 1; // Subtract 1 because element 0 is the function name + std::vector hyperparameters(num_args); + + for(int i = 0; i < num_args; ++i) { + SEXP arg = rhs[i + 1]; + if (Rcpp::is(arg) || Rcpp::is(arg)) { + // Evaluate the expression to get a numeric value + Rcpp::Environment base_env = Rcpp::Environment::base_env(); + SEXP result = Rcpp::Rcpp_eval(arg, base_env); + hyperparameters[i] = Rcpp::as(result); + } else { + hyperparameters[i] = Rcpp::as(arg); + } + } + + return FormulaComponents{ + transform, + fims::StringToDistributionLabel(Rcpp::as(distribution_string)), + hyperparameters + }; +} + +void setup_prior(Rcpp::Formula f, Rcpp::List parameter_vectors) { + + Rcpp::IntegerVector ids(parameter_vectors.size()); + + // Parse formula into components + FormulaComponents formula_parts = parse_distribution_formula(f); + fims::Distribution::Label distribution_name = formula_parts.distribution; + std::vector hyperparameters = formula_parts.hyperparameters; + fims::Transformation prior_transformation = formula_parts.transformation; + + // Pass prior transformation to ParameterVector + // Extract the raw pointer from the XPtr + for (int i = 0; i < parameter_vectors.size(); i++) { + ParameterVector pv = Rcpp::as(parameter_vectors[i]); + + // Modifies shared Transformation object + *pv.prior_transformation_m = prior_transformation; + ids[i] = pv.id_m; + } + + switch (distribution_name) { + + case fims::Distribution::Label::Normal: { + auto prior = std::make_shared(); + prior->expected_values[0].initial_value_m = hyperparameters[0]; + prior->log_sd[0].initial_value_m = fims_math::log(hyperparameters[1]); + prior->set_distribution_links("prior", ids); + break; + } + + case fims::Distribution::Label::Lognormal: { + auto prior = std::make_shared(); + prior->expected_values[0].initial_value_m = hyperparameters[0]; + prior->log_sd[0].initial_value_m = fims_math::log(hyperparameters[1]); + prior->set_distribution_links("prior", ids); + break; + } + + case fims::Distribution::Label::Gamma: { + auto prior = std::make_shared(); + prior->expected_values[0].initial_value_m = fims_math::log(hyperparameters[0]); + prior->log_sd[0].initial_value_m = fims_math::log(hyperparameters[1]); + prior->set_distribution_links("prior", ids); + break; + } + + case fims::Distribution::Label::InvGamma: { + auto prior = std::make_shared(); + prior->log_shape[0].initial_value_m = fims_math::log(hyperparameters[0]); + prior->log_scale[0].initial_value_m = fims_math::log(hyperparameters[1]); + prior->set_distribution_links("prior", ids); + break; + } + + default: + throw std::invalid_argument("Unsupported distribution type in add_prior."); + } +} + +//Define setup_prior_function in ParameterVector class +inline void register_prior_functions() { + ParameterVector::setup_prior_function = setup_prior; +} + +void add_shared_prior(Rcpp::Formula f, Rcpp::List parameter_vectors) { + setup_prior(f, parameter_vectors); +} + /** * @brief Clears the internal objects. * diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_distribution.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_distribution.hpp index 5de2996fa..e01bb19e4 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_distribution.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_distribution.hpp @@ -187,7 +187,7 @@ class DnormDistributionsInterface : public DistributionsInterfaceBase { */ ParameterVector expected_mean; /** - * @brief The uncertainty, which would be the standard deviation of x for the + * @brief The log_sd, which would be the standard deviation of x for the * normal distribution. */ ParameterVector log_sd; @@ -487,10 +487,10 @@ class DnormDistributionsInterface : public DistributionsInterfaceBase { info->RegisterParameter(distribution->log_sd[i]); } if (this->log_sd[i].estimation_type_m.get() == "random_effects") { - FIMS_ERROR_LOG("standard deviations cannot be set to random effects"); + FIMS_ERROR_LOG("log_sd values cannot be set to random effects"); } } - info->variable_map[this->log_sd.id_m] = &(distribution)->log_sd; + set_variable_map(&(distribution)->log_sd, this->log_sd); distribution->use_mean = this->use_mean_m.get(); distribution->expected_mean.resize(this->expected_mean.size()); @@ -507,8 +507,9 @@ class DnormDistributionsInterface : public DistributionsInterfaceBase { FIMS_ERROR_LOG("expected_mean cannot be set to random effects"); } } - info->variable_map[this->expected_mean.id_m] = - &(distribution)->expected_mean; + expected_mean.input_transformation_m->label = fims::Transformation::Label::identity; + expected_mean.prior_transformation_m->label = fims::Transformation::Label::identity; + set_variable_map(&(distribution)->expected_mean, this->expected_mean); info->density_components[distribution->id] = distribution; @@ -545,7 +546,7 @@ class DlnormDistributionsInterface : public DistributionsInterfaceBase { */ ParameterVector expected_values; /** - * @brief The uncertainty, which would be the natural logarithm of the + * @brief The log_sd, which would be the natural logarithm of the standard deviation (sd) of log(x) for this distribution. The natural log of the standard deviation is necessary because the exponential link function is applied to the log transformed standard deviation to insure @@ -829,7 +830,7 @@ class DlnormDistributionsInterface : public DistributionsInterfaceBase { FIMS_ERROR_LOG("standard deviations cannot be set to random effects"); } } - info->variable_map[this->log_sd.id_m] = &(distribution)->log_sd; + set_variable_map(&(distribution)->log_sd, this->log_sd); info->density_components[distribution->id] = distribution; diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_fleet.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_fleet.hpp index 2b508da05..376af3983 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_fleet.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_fleet.hpp @@ -440,7 +440,7 @@ class FleetInterface : public FleetInterfaceBase { } } // add to variable_map - info->variable_map[this->log_Fmort.id_m] = &(fleet)->log_Fmort; + info->variable_map[this->log_Fmort.id_m].variable = &(fleet)->log_Fmort; if (this->n_lengths.get() > 0) { fleet->age_to_length_conversion.resize( @@ -473,7 +473,7 @@ class FleetInterface : public FleetInterfaceBase { } } - info->variable_map[this->age_to_length_conversion.id_m] = + info->variable_map[this->age_to_length_conversion.id_m].variable = &(fleet)->age_to_length_conversion; } diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp index 8a83ab06f..f265e29bb 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp @@ -19,6 +19,7 @@ #include "../../../common/def.hpp" #include "../../../common/information.hpp" +#include "../../../common/types.hpp" #include "../../interface.hpp" #include "rcpp_shared_primitive.hpp" #include @@ -150,6 +151,16 @@ class ParameterVector { * @brief Parameter storage. */ std::shared_ptr> storage_m; + /** + * @brief The transformation of the parameter in the prior distribution. + */ + std::shared_ptr input_transformation_m = + std::make_shared(); + /** + * @brief The transformation of the parameter in the prior distribution. + */ + std::shared_ptr prior_transformation_m = + std::make_shared(); /** * @brief The local ID of the Parameter object. */ @@ -161,14 +172,19 @@ class ParameterVector { ParameterVector() { this->id_m = ParameterVector::id_g++; this->storage_m = std::make_shared>(); + this->input_transformation_m = std::make_shared(); + this->prior_transformation_m = std::make_shared(); this->storage_m->resize(1); // push_back(Rcpp::wrap(p)); } /** - * @brief The constructor. + * @brief The copy constructor. */ ParameterVector(const ParameterVector& other) - : storage_m(other.storage_m), id_m(other.id_m) {} + : storage_m(other.storage_m), + input_transformation_m(other.input_transformation_m), + prior_transformation_m(other.prior_transformation_m), + id_m(other.id_m) {} /** * @brief The constructor. @@ -176,6 +192,8 @@ class ParameterVector { ParameterVector(size_t size) { this->id_m = ParameterVector::id_g++; this->storage_m = std::make_shared>(); + this->input_transformation_m = std::make_shared(); + this->prior_transformation_m = std::make_shared(); this->storage_m->resize(size); for (size_t i = 0; i < size; i++) { storage_m->at(i) = Parameter(); @@ -197,6 +215,8 @@ class ParameterVector { this->storage_m = std::make_shared>(); // Use std::min to avoid comparing signed and unsigned types size_t n = std::min(static_cast(x.size()), size); + this->input_transformation_m = std::make_shared(); + this->prior_transformation_m = std::make_shared(); this->storage_m->resize(n); for (size_t i = 0; i < n; i++) { storage_m->at(i).initial_value_m = x[i]; @@ -211,6 +231,8 @@ class ParameterVector { ParameterVector(const fims::Vector& v) { this->id_m = ParameterVector::id_g++; this->storage_m = std::make_shared>(); + this->input_transformation_m = std::make_shared(); + this->prior_transformation_m = std::make_shared(); this->storage_m->resize(v.size()); for (size_t i = 0; i < v.size(); i++) { storage_m->at(i).initial_value_m = v[i]; @@ -250,6 +272,20 @@ class ParameterVector { return Rcpp::wrap(this->storage_m->at(pos - 1)); } + // register function first to avoid issues with circular dependancies + // in the include chain. The setup_prior_function will be set to equal the + // setup_prior() function in rcpp_interface.hpp. + static std::function setup_prior_function; + + void add_prior(Rcpp::Formula f) { + if (setup_prior_function) { + Rcpp::List pv_list = Rcpp::List::create(Rcpp::wrap(*this)); + setup_prior_function(f, pv_list); + } else { + throw std::runtime_error("setup_prior_function not registered."); + } + } + /** * @brief An internal accessor for calling a position of a ParameterVector * from R. @@ -346,7 +382,10 @@ class ParameterVector { } } }; -uint32_t ParameterVector::id_g = 0; +// static member definitions +uint32_t ParameterVector::id_g = 0; +std::function ParameterVector::setup_prior_function = nullptr; + /** * @brief Output for std::ostream& for a ParameterVector. @@ -400,7 +439,8 @@ class RealVector { * @brief The constructor. */ RealVector(const RealVector& other) - : storage_m(other.storage_m), id_m(other.id_m) {} + : storage_m(other.storage_m), + id_m(other.id_m) {} /** * @brief The constructor. @@ -645,8 +685,24 @@ class FIMSRcppInterfaceBase { } return ss.str(); } + +template +void set_variable_map(fims::Vector* ptr_fims_vector, + ParameterVector parameter_vector) { + + std::shared_ptr> info = + fims_info::Information::GetInstance(); + + info->variable_map[parameter_vector.id_m] = + typename fims_info::Information::VariableMapEntry( + ptr_fims_vector, + *parameter_vector.input_transformation_m, + *parameter_vector.prior_transformation_m); + } + }; std::vector> FIMSRcppInterfaceBase::fims_interface_objects; + #endif diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_models.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_models.hpp index c588e52a0..2d0680eb9 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_models.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_models.hpp @@ -1276,17 +1276,17 @@ class CatchAtAgeInterface : public FisheryModelInterfaceBase { fims::Vector{"n_years", "n_lengths"}); // replace elements in the variable map - info->variable_map[fleet_interface->log_landings_expected.id_m] = + info->variable_map[fleet_interface->log_landings_expected.id_m].variable = &(derived_quantities["log_landings_expected"]); - info->variable_map[fleet_interface->log_index_expected.id_m] = + info->variable_map[fleet_interface->log_index_expected.id_m].variable = &(derived_quantities["log_index_expected"]); - info->variable_map[fleet_interface->agecomp_expected.id_m] = + info->variable_map[fleet_interface->agecomp_expected.id_m].variable = &(derived_quantities["agecomp_expected"]); - info->variable_map[fleet_interface->agecomp_proportion.id_m] = + info->variable_map[fleet_interface->agecomp_proportion.id_m].variable = &(derived_quantities["agecomp_proportion"]); - info->variable_map[fleet_interface->lengthcomp_expected.id_m] = + info->variable_map[fleet_interface->lengthcomp_expected.id_m].variable = &(derived_quantities["lengthcomp_expected"]); - info->variable_map[fleet_interface->lengthcomp_proportion.id_m] = + info->variable_map[fleet_interface->lengthcomp_proportion.id_m].variable = &(derived_quantities["lengthcomp_proportion"]); } diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp index edb73f1af..7ed5806dc 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp @@ -356,7 +356,7 @@ class PopulationInterface : public PopulationInterfaceBase { population->spawning_biomass_ratio.resize( this->spawning_biomass_ratio.size()); } - info->variable_map[this->spawning_biomass_ratio.id_m] = + info->variable_map[this->spawning_biomass_ratio.id_m].variable = &(population)->spawning_biomass_ratio; population->log_init_naa.resize(this->log_init_naa.size()); @@ -375,7 +375,7 @@ class PopulationInterface : public PopulationInterfaceBase { info->RegisterRandomEffect(population->log_M[i]); } } - info->variable_map[this->log_M.id_m] = &(population)->log_M; + info->variable_map[this->log_M.id_m].variable = &(population)->log_M; for (size_t i = 0; i < log_f_multiplier.size(); i++) { population->log_f_multiplier[i] = @@ -397,7 +397,7 @@ class PopulationInterface : public PopulationInterfaceBase { info->RegisterRandomEffect(population->log_f_multiplier[i]); } } - info->variable_map[this->log_f_multiplier.id_m] = + info->variable_map[this->log_f_multiplier.id_m].variable = &(population)->log_f_multiplier; for (size_t i = 0; i < log_init_naa.size(); i++) { @@ -417,7 +417,7 @@ class PopulationInterface : public PopulationInterfaceBase { info->RegisterRandomEffect(population->log_init_naa[i]); } } - info->variable_map[this->log_init_naa.id_m] = &(population)->log_init_naa; + info->variable_map[this->log_init_naa.id_m].variable = &(population)->log_init_naa; for (size_t i = 0; i < ages.size(); i++) { population->ages[i] = this->ages[i]; diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp index e2ed1ccd9..9251b533c 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp @@ -353,7 +353,14 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { info->RegisterRandomEffect(recruitment->logit_steep[i]); } } - info->variable_map[this->logit_steep.id_m] = &(recruitment)->logit_steep; + + //set transformations for parameter since default is log + logit_steep.input_transformation_m->label = fims::Transformation::Label::logit; + logit_steep.input_transformation_m->args.lower = 0.2; + logit_steep.input_transformation_m->args.upper = 1.0; + logit_steep.prior_transformation_m->label = fims::Transformation::Label::logit; + logit_steep.prior_transformation_m->args.lower = 0.2; + logit_steep.prior_transformation_m->args.upper = 1.0; // set log_rzero recruitment->log_rzero.resize(this->log_rzero.size()); @@ -375,7 +382,7 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { info->RegisterRandomEffect(recruitment->log_rzero[i]); } } - info->variable_map[this->log_rzero.id_m] = &(recruitment)->log_rzero; + info->variable_map[this->log_rzero.id_m].variable = &(recruitment)->log_rzero; // set log_recruit_devs recruitment->log_recruit_devs.resize(this->log_devs.size()); for (size_t i = 0; i < this->log_devs.size(); i++) { @@ -397,7 +404,7 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { } } - info->variable_map[this->log_devs.id_m] = &(recruitment)->log_recruit_devs; + info->variable_map[this->log_devs.id_m].variable = &(recruitment)->log_recruit_devs; // set log_r recruitment->log_r.resize(this->log_r.size()); @@ -418,13 +425,13 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { } } - info->variable_map[this->log_r.id_m] = &(recruitment)->log_r; + info->variable_map[this->log_r.id_m].variable = &(recruitment)->log_r; // set log_expected_recruitment recruitment->log_expected_recruitment.resize(this->n_years.get() - 1); for (size_t i = 0; i < static_cast(this->n_years.get() - 1); i++) { recruitment->log_expected_recruitment[i] = 0; } - info->variable_map[this->log_expected_recruitment.id_m] = + info->variable_map[this->log_expected_recruitment.id_m].variable = &(recruitment)->log_expected_recruitment; // add to Information diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp index 5656940ea..ab9c675b7 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp @@ -255,7 +255,7 @@ class LogisticSelectivityInterface : public SelectivityInterfaceBase { info->RegisterRandomEffectName(ss.str()); } } - info->variable_map[this->inflection_point.id_m] = + info->variable_map[this->inflection_point.id_m].variable = &(selectivity)->inflection_point; selectivity->slope.resize(this->slope.size()); @@ -274,7 +274,7 @@ class LogisticSelectivityInterface : public SelectivityInterfaceBase { info->RegisterRandomEffect(selectivity->slope[i]); } } - info->variable_map[this->slope.id_m] = &(selectivity)->slope; + info->variable_map[this->slope.id_m].variable = &(selectivity)->slope; // add to Information info->selectivity_models[selectivity->id] = selectivity; @@ -514,7 +514,7 @@ class DoubleLogisticSelectivityInterface : public SelectivityInterfaceBase { info->RegisterRandomEffect(selectivity->inflection_point_asc[i]); } } - info->variable_map[this->inflection_point_asc.id_m] = + info->variable_map[this->inflection_point_asc.id_m].variable = &(selectivity)->inflection_point_asc; selectivity->slope_asc.resize(this->slope_asc.size()); @@ -536,7 +536,7 @@ class DoubleLogisticSelectivityInterface : public SelectivityInterfaceBase { info->RegisterRandomEffect(selectivity->slope_asc[i]); } } - info->variable_map[this->slope_asc.id_m] = &(selectivity)->slope_asc; + info->variable_map[this->slope_asc.id_m].variable = &(selectivity)->slope_asc; selectivity->inflection_point_desc.resize( this->inflection_point_desc.size()); @@ -561,7 +561,7 @@ class DoubleLogisticSelectivityInterface : public SelectivityInterfaceBase { info->RegisterRandomEffect(selectivity->inflection_point_desc[i]); } } - info->variable_map[this->inflection_point_desc.id_m] = + info->variable_map[this->inflection_point_desc.id_m].variable = &(selectivity)->inflection_point_desc; selectivity->slope_desc.resize(this->slope_desc.size()); @@ -584,7 +584,7 @@ class DoubleLogisticSelectivityInterface : public SelectivityInterfaceBase { } } - info->variable_map[this->slope_desc.id_m] = &(selectivity)->slope_desc; + info->variable_map[this->slope_desc.id_m].variable = &(selectivity)->slope_desc; // add to Information info->selectivity_models[selectivity->id] = selectivity; diff --git a/man/Cpp_functions.Rd b/man/Cpp_functions.Rd index d83bfecbe..ad8b0c0a4 100644 --- a/man/Cpp_functions.Rd +++ b/man/Cpp_functions.Rd @@ -28,6 +28,7 @@ C++ docs. \details{ \itemize{ \item \href{https://noaa-fims.github.io/FIMS/doxygen/rcpp__interface_8hpp.html}{clear} +\item \href{https://noaa-fims.github.io/FIMS/doxygen/add__shared__prior_8hpp.html}{add_shared_prior} \item \href{https://noaa-fims.github.io/FIMS/doxygen/rcpp__interface_8hpp.html}{get_fixed} \item \href{https://noaa-fims.github.io/FIMS/doxygen/rcpp__interface_8hpp.html}{get_log} \item \href{https://noaa-fims.github.io/FIMS/doxygen/rcpp__interface_8hpp.html}{get_log_errors} diff --git a/src/FIMS.cpp b/src/FIMS.cpp index 6791ba5df..61f1792c3 100644 --- a/src/FIMS.cpp +++ b/src/FIMS.cpp @@ -16,6 +16,8 @@ template Type objective_function::operator()() { + DATA_INTEGER(do_mcmc); + PARAMETER_VECTOR(p); PARAMETER_VECTOR(re); @@ -38,6 +40,7 @@ Type objective_function::operator()() { *information->random_effects_parameters[i] = re[i]; } model -> of = this; + model -> jacobian_flag = (do_mcmc == 1); Type nll = 0; //evaluate the model objective function value diff --git a/src/fims_modules.hpp b/src/fims_modules.hpp index 77987dace..6e25c9225 100644 --- a/src/fims_modules.hpp +++ b/src/fims_modules.hpp @@ -89,6 +89,9 @@ RCPP_MODULE(fims) { "get_random_names", &get_random_names, "See " "https://noaa-fims.github.io/FIMS/doxygen/rcpp__interface_8hpp.html."); + Rcpp::function("add_shared_prior", &setup_prior, + "See " + "https://noaa-fims.github.io/FIMS/doxygen/rcpp__interface_8hpp.html."); Rcpp::function( "clear", clear, "See " @@ -165,6 +168,7 @@ RCPP_MODULE(fims) { .constructor() .constructor() .constructor() + .method("add_prior", &ParameterVector::add_prior) .method("get", &ParameterVector::get) .method("set", &ParameterVector::set) .method("show", &ParameterVector::show) diff --git a/src/init.hpp b/src/init.hpp index 905da2b18..1e28ebe10 100644 --- a/src/init.hpp +++ b/src/init.hpp @@ -35,6 +35,7 @@ static const R_CallMethodDef CallEntries[] = { void R_init_FIMS(DllInfo *dll) { R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); R_useDynamicSymbols(dll, FALSE); + register_prior_functions(); #ifdef TMB_CCALLABLES TMB_CCALLABLES("FIMS"); #endif diff --git a/tests/gtest/CMakeLists.txt b/tests/gtest/CMakeLists.txt index 354e58bfe..ff9b8b00d 100644 --- a/tests/gtest/CMakeLists.txt +++ b/tests/gtest/CMakeLists.txt @@ -278,6 +278,30 @@ target_link_libraries(fims_math_logit gtest_discover_tests(fims_math_logit) +# test_fims_transformation_ApplyTransformation.cpp +add_executable(fims_apply_transformation + test_fims_transformations_ApplyTransformation.cpp +) + +target_link_libraries(fims_apply_transformation + gtest_main + fims_test +) + +gtest_discover_tests(fims_apply_transformation) + +# test_fims_transformation_TransformPrior.cpp +add_executable(fims_transform_prior + test_fims_transformations_TransformPrior.cpp +) + +target_link_libraries(fims_transform_prior + gtest_main + fims_test +) + +gtest_discover_tests(fims_transform_prior) + # test_DoubleLogistic_DoubleLogisticSelectivity_Evaluate.cpp add_executable(DoubleLogistic_DoubleLogisticSelectivity_Evaluate test_DoubleLogistic_DoubleLogisticSelectivity_Evaluate.cpp diff --git a/tests/gtest/test_fims_transformations_ApplyTransformation.cpp b/tests/gtest/test_fims_transformations_ApplyTransformation.cpp new file mode 100644 index 000000000..cc5825179 --- /dev/null +++ b/tests/gtest/test_fims_transformations_ApplyTransformation.cpp @@ -0,0 +1,244 @@ +#include "gtest/gtest.h" +#include "common/fims_transformations.hpp" + +namespace { + + +// ============================================================ +// ApplyTransformation Tests +// ============================================================ + + +TEST(ApplyTransformation, IdentityReturnsInputUnchanged) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::identity; + + EXPECT_EQ(fims_transformations::ApplyTransformation(0.0, trans), 0.0); + EXPECT_EQ(fims_transformations::ApplyTransformation(1.0, trans), 1.0); + EXPECT_EQ(fims_transformations::ApplyTransformation(-1.0, trans), -1.0); + EXPECT_EQ(fims_transformations::ApplyTransformation(100.0, trans), 100.0); +} + +TEST(ApplyTransformation, LogTransformMatchesStdLog) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::log; + + EXPECT_EQ(fims_transformations::ApplyTransformation(1.0, trans), std::log(1.0)); + EXPECT_EQ(fims_transformations::ApplyTransformation(10.0, trans), std::log(10.0)); + EXPECT_EQ(fims_transformations::ApplyTransformation(0.5, trans), std::log(0.5)); + EXPECT_EQ(fims_transformations::ApplyTransformation(100.0, trans), std::log(100.0)); +} + +TEST(ApplyTransformation, ExpTransformMatchesStdExp) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::exp; + + EXPECT_EQ(fims_transformations::ApplyTransformation(0.0, trans), std::exp(0.0)); + EXPECT_EQ(fims_transformations::ApplyTransformation(1.0, trans), std::exp(1.0)); + EXPECT_EQ(fims_transformations::ApplyTransformation(-1.0, trans), std::exp(-1.0)); + EXPECT_EQ(fims_transformations::ApplyTransformation(3.0, trans), std::exp(3.0)); +} + +TEST(ApplyTransformation, SquareTransformReturnsSquaredValue) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::square; + + EXPECT_EQ(fims_transformations::ApplyTransformation(2.0, trans), 4.0); + EXPECT_EQ(fims_transformations::ApplyTransformation(3.0, trans), 9.0); + EXPECT_EQ(fims_transformations::ApplyTransformation(0.5, trans), 0.25); + EXPECT_EQ(fims_transformations::ApplyTransformation(0.0, trans), 0.0); +} + +TEST(ApplyTransformation, SqrtTransformMatchesStdSqrt) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::sqrt; + + EXPECT_EQ(fims_transformations::ApplyTransformation(4.0, trans), std::sqrt(4.0)); + EXPECT_EQ(fims_transformations::ApplyTransformation(9.0, trans), std::sqrt(9.0)); + EXPECT_EQ(fims_transformations::ApplyTransformation(0.25, trans), std::sqrt(0.25)); + EXPECT_EQ(fims_transformations::ApplyTransformation(0.0, trans), 0.0); +} + +TEST(ApplyTransformation, LogitTransformDefaultBounds) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::logit; + // Default bounds: lower = 0.0, upper = 1.0 + + // logit(0.5) = 0 + EXPECT_EQ(fims_transformations::ApplyTransformation(0.5, trans), 0.0); + // logit(0.731) ~ 1.0 + EXPECT_NEAR(fims_transformations::ApplyTransformation(0.731059, trans), 1.0, 1e-5); + // logit(0.269) ~ -1.0 + EXPECT_NEAR(fims_transformations::ApplyTransformation(0.268941, trans), -1.0, 1e-5); +} + +TEST(ApplyTransformation, LogitTransformCustomBounds) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::logit; + trans.args.lower = 0.2; + trans.args.upper = 1.0; + + // At midpoint of [0.2, 1.0] = 0.6, logit should be ~ 0 + EXPECT_NEAR(fims_transformations::ApplyTransformation(0.6, trans), 0.0, 1e-5); +} + +TEST(ApplyTransformation, UnsupportedLabelThrows) { + fims::Transformation trans; + // Force an invalid label value + trans.label = static_cast(999); + + EXPECT_THROW( + fims_transformations::ApplyTransformation(1.0, trans), + std::invalid_argument + ); +} + + +// ============================================================ +// ApplyBackTransformation Tests +// ============================================================ + + + +TEST(ApplyBackTransformation, IdentityReturnsInputUnchanged) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::identity; + + EXPECT_EQ(fims_transformations::ApplyBackTransformation(0.0, trans), 0.0); + EXPECT_EQ(fims_transformations::ApplyBackTransformation(1.0, trans), 1.0); + EXPECT_EQ(fims_transformations::ApplyBackTransformation(-1.0, trans), -1.0); + EXPECT_EQ(fims_transformations::ApplyBackTransformation(100.0, trans), 100.0); +} + +TEST(ApplyBackTransformation, LogBackTransformIsExp) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::log; + + // Back transform of log is exp + EXPECT_EQ(fims_transformations::ApplyBackTransformation(0.0, trans), std::exp(0.0)); + EXPECT_EQ(fims_transformations::ApplyBackTransformation(1.0, trans), std::exp(1.0)); + EXPECT_EQ(fims_transformations::ApplyBackTransformation(-1.0, trans), std::exp(-1.0)); + EXPECT_EQ(fims_transformations::ApplyBackTransformation(3.0, trans), std::exp(3.0)); +} + +TEST(ApplyBackTransformation, ExpBackTransformIsLog) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::exp; + + // Back transform of exp is log + EXPECT_EQ(fims_transformations::ApplyBackTransformation(1.0, trans), std::log(1.0)); + EXPECT_EQ(fims_transformations::ApplyBackTransformation(10.0, trans), std::log(10.0)); + EXPECT_EQ(fims_transformations::ApplyBackTransformation(0.5, trans), std::log(0.5)); +} + +TEST(ApplyBackTransformation, SquareBackTransformIsSqrt) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::square; + + // Back transform of square is sqrt + EXPECT_EQ(fims_transformations::ApplyBackTransformation(4.0, trans), 2.0); + EXPECT_EQ(fims_transformations::ApplyBackTransformation(9.0, trans), 3.0); + EXPECT_EQ(fims_transformations::ApplyBackTransformation(0.25, trans), 0.5); +} + +TEST(ApplyBackTransformation, SqrtBackTransformIsSquare) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::sqrt; + + // Back transform of sqrt is square + EXPECT_EQ(fims_transformations::ApplyBackTransformation(2.0, trans), 4.0); + EXPECT_EQ(fims_transformations::ApplyBackTransformation(3.0, trans), 9.0); + EXPECT_EQ(fims_transformations::ApplyBackTransformation(0.5, trans), 0.25); +} + +TEST(ApplyBackTransformation, LogitBackTransformIsInvLogit) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::logit; + // Default bounds: lower = 0.0, upper = 1.0 + + // inv_logit(0) = 0.5 + EXPECT_EQ(fims_transformations::ApplyBackTransformation(0.0, trans), 0.5); + // inv_logit(1) ~ 0.731 + EXPECT_NEAR(fims_transformations::ApplyBackTransformation(1.0, trans), 0.731059, 1e-5); + // inv_logit(-1) ~ 0.269 + EXPECT_NEAR(fims_transformations::ApplyBackTransformation(-1.0, trans), 0.268941, 1e-5); +} + +TEST(ApplyBackTransformation, UnsupportedLabelThrows) { + fims::Transformation trans; + trans.label = static_cast(999); + + EXPECT_THROW( + fims_transformations::ApplyBackTransformation(1.0, trans), + std::invalid_argument + ); +} + + +// ============================================================ +// Round-trip tests: ApplyTransformation(ApplyBackTransformation(x)) == x +// ============================================================ + +TEST(TransformationRoundTrip, LogAndBackLog) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::log; + + double x = 5.0; + double forward = fims_transformations::ApplyTransformation(x, trans); + double back = fims_transformations::ApplyBackTransformation(forward, trans); + EXPECT_NEAR(back, x, 1e-8); +} + +TEST(TransformationRoundTrip, ExpAndBackExp) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::exp; + + double x = 2.0; + double forward = fims_transformations::ApplyTransformation(x, trans); + double back = fims_transformations::ApplyBackTransformation(forward, trans); + EXPECT_NEAR(back, x, 1e-8); +} + +TEST(TransformationRoundTrip, SquareAndBackSquare) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::square; + + double x = 3.0; + double forward = fims_transformations::ApplyTransformation(x, trans); + double back = fims_transformations::ApplyBackTransformation(forward, trans); + EXPECT_NEAR(back, x, 1e-8); +} + +TEST(TransformationRoundTrip, SqrtAndBackSqrt) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::sqrt; + + double x = 7.0; + double forward = fims_transformations::ApplyTransformation(x, trans); + double back = fims_transformations::ApplyBackTransformation(forward, trans); + EXPECT_NEAR(back, x, 1e-8); +} + +TEST(TransformationRoundTrip, LogitAndBackLogitDefaultBounds) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::logit; + + double x = 0.7; + double forward = fims_transformations::ApplyTransformation(x, trans); + double back = fims_transformations::ApplyBackTransformation(forward, trans); + EXPECT_NEAR(back, x, 1e-8); +} + +TEST(TransformationRoundTrip, LogitAndBackLogitCustomBounds) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::logit; + trans.args.lower = 0.2; + trans.args.upper = 1.0; + + double x = 0.8; + double forward = fims_transformations::ApplyTransformation(x, trans); + double back = fims_transformations::ApplyBackTransformation(forward, trans); + EXPECT_NEAR(back, x, 1e-8); +} + +} \ No newline at end of file diff --git a/tests/gtest/test_fims_transformations_TransformPrior.cpp b/tests/gtest/test_fims_transformations_TransformPrior.cpp new file mode 100644 index 000000000..47b516734 --- /dev/null +++ b/tests/gtest/test_fims_transformations_TransformPrior.cpp @@ -0,0 +1,134 @@ +#include "gtest/gtest.h" +#include "common/fims_transformations.hpp" + +namespace { + + +// ============================================================ +// TransformPrior (vector overload) Tests +// ============================================================ + +TEST(TransformPriorVector, IdenticalTransformationsReturnInputUnchanged) { + fims::Transformation trans; + trans.label = fims::Transformation::Label::log; + + fims::Vector input; + input.resize(3); + input[0] = 1.0; + input[1] = 2.0; + input[2] = 3.0; + + // Same input and prior transformation - should return input unchanged + fims::Vector result = fims_transformations::TransformPrior(input, trans, trans); + EXPECT_NEAR(result[0], 1.0, 1e-8); + EXPECT_NEAR(result[1], 2.0, 1e-8); + EXPECT_NEAR(result[2], 3.0, 1e-8); +} + +TEST(TransformPriorVector, LogInputToNaturalScale) { + fims::Transformation log_trans; + log_trans.label = fims::Transformation::Label::log; + + fims::Transformation identity_trans; + identity_trans.label = fims::Transformation::Label::identity; + + // Input is log(sd), prior is on natural scale (sd) + fims::Vector input; + input.resize(3); + input[0] = std::log(1.0); + input[1] = std::log(2.0); + input[2] = std::log(3.0); + + fims::Vector result = fims_transformations::TransformPrior( + input, log_trans, identity_trans); + + EXPECT_NEAR(result[0], 1.0, 1e-8); + EXPECT_NEAR(result[1], 2.0, 1e-8); + EXPECT_NEAR(result[2], 3.0, 1e-8); +} + +TEST(TransformPriorVector, LogInputToVarianceScale) { + fims::Transformation log_trans; + log_trans.label = fims::Transformation::Label::log; + + fims::Transformation square_trans; + square_trans.label = fims::Transformation::Label::square; + + // Input is log(sd), prior is on variance (sd^2) + // log(sd) -> exp -> sd -> square -> sd^2 + double sd = 2.0; + fims::Vector input; + input.resize(1); + input[0] = std::log(sd); + + fims::Vector result = fims_transformations::TransformPrior( + input, log_trans, square_trans); + + EXPECT_NEAR(result[0], sd * sd, 1e-8); +} + +TEST(TransformPriorVector, PreservesVectorSize) { + fims::Transformation log_trans; + log_trans.label = fims::Transformation::Label::log; + + fims::Transformation identity_trans; + identity_trans.label = fims::Transformation::Label::identity; + + fims::Vector input; + input.resize(5); + for (size_t i = 0; i < 5; i++) input[i] = std::log(i + 1.0); + + fims::Vector result = fims_transformations::TransformPrior( + input, log_trans, identity_trans); + + EXPECT_EQ(5, result.size()); +} + +// ============================================================ +// TransformPrior (scalar overload) Tests +// ============================================================ + +TEST(TransformPriorScalar, LogInputToNaturalScale) { + fims::Transformation log_trans; + log_trans.label = fims::Transformation::Label::log; + + fims::Transformation identity_trans; + identity_trans.label = fims::Transformation::Label::identity; + + double log_sd = std::log(3.0); + double result = fims_transformations::TransformPrior( + log_sd, log_trans, identity_trans); + + EXPECT_NEAR(result, 3.0, 1e-8); +} + +TEST(TransformPriorScalar, LogInputToVarianceScale) { + fims::Transformation log_trans; + log_trans.label = fims::Transformation::Label::log; + + fims::Transformation square_trans; + square_trans.label = fims::Transformation::Label::square; + + double sd = 2.5; + double result = fims_transformations::TransformPrior( + std::log(sd), log_trans, square_trans); + + EXPECT_NEAR(result, sd * sd, 1e-8); +} + +TEST(TransformPriorScalar, IdentityInputToSquareScale) { + fims::Transformation identity_trans; + identity_trans.label = fims::Transformation::Label::identity; + + fims::Transformation square_trans; + square_trans.label = fims::Transformation::Label::square; + + double x = 4.0; + double result = fims_transformations::TransformPrior( + x, identity_trans, square_trans); + + EXPECT_NEAR(result, 16.0, 1e-8); +} + + +} \ No newline at end of file diff --git a/tests/gtest/test_info_setup_priors.cpp b/tests/gtest/test_info_setup_priors.cpp index d213e4b09..124f40f83 100644 --- a/tests/gtest/test_info_setup_priors.cpp +++ b/tests/gtest/test_info_setup_priors.cpp @@ -27,10 +27,10 @@ namespace selectivity2->slope[0] = 0.18; // Set up variable map to point to selectivity parameters - info->variable_map[1] = &(selectivity1)->inflection_point; - info->variable_map[2] = &(selectivity1)->slope; - info->variable_map[3] = &(selectivity2)->inflection_point; - info->variable_map[4] = &(selectivity2)->slope; + info->variable_map[1].variable = &(selectivity1)->inflection_point; + info->variable_map[2].variable = &(selectivity1)->slope; + info->variable_map[3].variable = &(selectivity2)->inflection_point; + info->variable_map[4].variable = &(selectivity2)->slope; //Create new normal distributions std::shared_ptr > normal_inflection_point = @@ -54,6 +54,8 @@ namespace // Call function that links key ID to variable map pointers given variable map ID // This function will set the density component, priors, to the respective parameters info->SetupPriors(); + normal_inflection_point->Prepare(); + lognormal_slope->Prepare(); EXPECT_EQ((*normal_inflection_point->priors[0])[0], selectivity1->inflection_point[0]); EXPECT_EQ((*normal_inflection_point->priors[1])[0], selectivity2->inflection_point[0]); @@ -61,8 +63,8 @@ namespace EXPECT_EQ((*lognormal_slope->priors[1])[0], selectivity2->slope[0]); EXPECT_EQ((normal_inflection_point->get_observed(0)), selectivity1->inflection_point[0]); EXPECT_EQ((normal_inflection_point->get_observed(1)), selectivity2->inflection_point[0]); - EXPECT_EQ((lognormal_slope->get_observed(0)), selectivity1->slope[0]); - EXPECT_EQ((lognormal_slope->get_observed(1)), selectivity2->slope[0]); + EXPECT_NEAR((lognormal_slope->get_observed(0)), selectivity1->slope[0], 1e-10); + EXPECT_NEAR((lognormal_slope->get_observed(1)), selectivity2->slope[0], 1e-10); //update the value in normal to check if it is updated in selectivity (*normal_inflection_point->priors[0])[0] = 20.5; diff --git a/tests/gtest/test_info_setup_random_effects.cpp b/tests/gtest/test_info_setup_random_effects.cpp index 99b740f24..351067d68 100644 --- a/tests/gtest/test_info_setup_random_effects.cpp +++ b/tests/gtest/test_info_setup_random_effects.cpp @@ -19,8 +19,8 @@ namespace recruitment->log_expected_recruitment.resize(2); recruitment->log_expected_recruitment[0] = 2.1; recruitment->log_expected_recruitment[1] = -1.7; - info->variable_map[1] = &(recruitment)->log_r; - info->variable_map[2] = &(recruitment)->log_expected_recruitment; + info->variable_map[1].variable = &(recruitment)->log_r; + info->variable_map[2].variable = &(recruitment)->log_expected_recruitment; //Create a new normal distribution std::shared_ptr > normal = diff --git a/tests/testthat/helper-integration-tests-setup-function.R b/tests/testthat/helper-integration-tests-setup-function.R index 158be443d..3b2a23dd8 100644 --- a/tests/testthat/helper-integration-tests-setup-function.R +++ b/tests/testthat/helper-integration-tests-setup-function.R @@ -367,6 +367,7 @@ setup_and_run_FIMS_without_wrappers <- function(iter_id, recruitment_distribution$log_sd$resize(1) recruitment_distribution$log_sd[1]$value <- log(om_input[["logR_sd"]]) + recruitment_distribution$observed_values$resize(om_input[["nyr"]] - 1) recruitment_distribution$expected_values$resize(om_input[["nyr"]] - 1) for (i in 1:(om_input[["nyr"]] - 1)) { diff --git a/tests/testthat/test-integration-fims-bayesian-prior-predictive.R b/tests/testthat/test-integration-fims-bayesian-prior-predictive.R index 6df7f6bda..074c5d79a 100644 --- a/tests/testthat/test-integration-fims-bayesian-prior-predictive.R +++ b/tests/testthat/test-integration-fims-bayesian-prior-predictive.R @@ -148,23 +148,23 @@ test_that("posterior equals prior with no data", { # Set up priors for selectivity parameters and link to both fishery and survey selectivity slope_mean <- mean(c(om_input[["sel_fleet"]][["fleet1"]][["slope.sel1"]], om_input[["sel_survey"]][["survey1"]][["slope.sel1"]])) - slope_sd <- 3 + slope_var <- 9 slope_prior <- methods::new(DnormDistribution) slope_prior$expected_values$resize(2) slope_prior$expected_values[1]$value <- slope_mean slope_prior$expected_values[2]$value <- slope_mean slope_prior$log_sd$resize(1) - slope_prior$log_sd[1]$value <- log(slope_sd) + slope_prior$log_sd[1]$value <- slope_var slope_prior$set_distribution_links("prior", c(fishing_fleet_selectivity$slope$get_id(), survey_fleet_selectivity$slope$get_id())) inflection_point_mean <- mean(c(om_input[["sel_fleet"]][["fleet1"]][["A50.sel1"]], om_input[["sel_survey"]][["survey1"]][["A50.sel1"]])) - inflection_point_sd <- 3 + inflection_point_var <- 9 inflection_point_prior <- methods::new(DnormDistribution) inflection_point_prior$expected_values$resize(2) inflection_point_prior$expected_values[1]$value <- inflection_point_mean inflection_point_prior$expected_values[2]$value <- inflection_point_mean inflection_point_prior$log_sd$resize(1) - inflection_point_prior$log_sd[1]$value <- log(inflection_point_sd) + inflection_point_prior$log_sd[1]$value <- inflection_point_var inflection_point_prior$set_distribution_links("prior", c(fishing_fleet_selectivity$inflection_point$get_id(), survey_fleet_selectivity$inflection_point$get_id())) @@ -251,11 +251,11 @@ test_that("posterior equals prior with no data", { slope_input <- c(om_input[["sel_fleet"]][["fleet1"]][["slope.sel1"]], om_input[["sel_survey"]][["survey1"]][["slope.sel1"]]) #' @description Test the slope nll expect_equal( - report_nll[1], -sum(dnorm(slope_input, mean = slope_mean, sd = 3, log = TRUE)) + report_nll[1], -sum(dnorm(slope_input, mean = slope_mean, sd = sqrt(slope_var), log = TRUE)) ) #' @description Test the inflection point nll expect_equal( - report_nll[2], -sum(dnorm(inflection_point_input, mean = inflection_point_mean, sd = 3, log = TRUE)) + report_nll[2], -sum(dnorm(inflection_point_input, mean = inflection_point_mean, sd = sqrt(inflection_point_var), log = TRUE)) ) # Fit MCMC using SparseNUTS @@ -282,9 +282,9 @@ test_that("posterior equals prior with no data", { #' @description Test that the posterior means for slope match the prior means. expect_equal(slope_est[[i]], slope_mean) #' @description Test that the posterior standard errors for inflection point match the prior standard errors. - expect_equal(inflection_point_se[[i]], inflection_point_sd) + expect_equal(inflection_point_se[[i]], sqrt(inflection_point_var)) #' @description Test that the posterior standard errors for slope match the prior standard errors. - expect_equal(slope_se[[i]], slope_sd) + expect_equal(slope_se[[i]], sqrt(slope_var)) } clear() diff --git a/tests/testthat/test-rcpp-distribution.R b/tests/testthat/test-rcpp-distribution.R index 952ba46c3..2ce85dce8 100644 --- a/tests/testthat/test-rcpp-distribution.R +++ b/tests/testthat/test-rcpp-distribution.R @@ -250,6 +250,7 @@ test_that("rcpp_distribution returns correct outputs for edge cases", { dnorm_$observed_values[1]$value <- y dnorm_$expected_values[1]$value <- 0 dnorm_$log_sd[1]$value <- -10 + # evaluate the density and compare with R #' @description Test extreme log_sd values for dnorm (-10) return expected output. expect_equal(dnorm_$evaluate(), stats::dnorm(y, 0, exp(-10), TRUE)) From e8a084c0684e233c0f9430158993edea68ca8391 Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA Date: Fri, 29 May 2026 17:03:11 -0700 Subject: [PATCH 2/4] remove gamma and invgamma from distribution enum --- inst/include/interface/rcpp/rcpp_interface.hpp | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/inst/include/interface/rcpp/rcpp_interface.hpp b/inst/include/interface/rcpp/rcpp_interface.hpp index 993fbae69..4863703b2 100644 --- a/inst/include/interface/rcpp/rcpp_interface.hpp +++ b/inst/include/interface/rcpp/rcpp_interface.hpp @@ -421,22 +421,6 @@ void setup_prior(Rcpp::Formula f, Rcpp::List parameter_vectors) { break; } - case fims::Distribution::Label::Gamma: { - auto prior = std::make_shared(); - prior->expected_values[0].initial_value_m = fims_math::log(hyperparameters[0]); - prior->log_sd[0].initial_value_m = fims_math::log(hyperparameters[1]); - prior->set_distribution_links("prior", ids); - break; - } - - case fims::Distribution::Label::InvGamma: { - auto prior = std::make_shared(); - prior->log_shape[0].initial_value_m = fims_math::log(hyperparameters[0]); - prior->log_scale[0].initial_value_m = fims_math::log(hyperparameters[1]); - prior->set_distribution_links("prior", ids); - break; - } - default: throw std::invalid_argument("Unsupported distribution type in add_prior."); } From fb03fa296df53a3e6a190b10b35617f2c5a8e55a Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA Date: Fri, 29 May 2026 17:21:55 -0700 Subject: [PATCH 3/4] cleanup - remove changes that were pulled over from surplus production --- CMakeLists.txt | 9 +-------- inst/include/common/fims_vector.hpp | 1 - inst/include/distributions/functors/normal_lpdf.hpp | 2 +- .../interface/rcpp/rcpp_objects/rcpp_distribution.hpp | 6 +++--- .../interface/rcpp/rcpp_objects/rcpp_interface_base.hpp | 4 +--- tests/testthat/helper-integration-tests-setup-function.R | 1 - tests/testthat/test-rcpp-distribution.R | 1 - 7 files changed, 6 insertions(+), 18 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2a8a17079..c6d1cbe1d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -121,19 +121,12 @@ enable_testing() # Creates an interface library fims_test add_library(fims_test INTERFACE) -# Get the exact include path directly from R -EXECUTE_PROCESS( - COMMAND R "--slave" "--no-save" "-e" "cat(R.home('include'))" - OUTPUT_VARIABLE R_INCLUDE_DIR - OUTPUT_STRIP_TRAILING_WHITESPACE -) - # Add include directories to the fims_test target target_include_directories(fims_test INTERFACE inst/include ${rcpp_SOURCE_DIR}/inst/include - ${R_INCLUDE_DIR} + ${R_HOME}/include ) # Add compile definition STD_LIB to the fims_test target. diff --git a/inst/include/common/fims_vector.hpp b/inst/include/common/fims_vector.hpp index 6b899afd3..176ae7c11 100644 --- a/inst/include/common/fims_vector.hpp +++ b/inst/include/common/fims_vector.hpp @@ -8,7 +8,6 @@ #ifndef FIMS_VECTOR_HPP #define FIMS_VECTOR_HPP -#include "types.hpp" #include "../interface/interface.hpp" #include #include diff --git a/inst/include/distributions/functors/normal_lpdf.hpp b/inst/include/distributions/functors/normal_lpdf.hpp index 24238ddfa..a0b60aa91 100644 --- a/inst/include/distributions/functors/normal_lpdf.hpp +++ b/inst/include/distributions/functors/normal_lpdf.hpp @@ -114,7 +114,7 @@ struct NormalLPDF : public DensityComponentBase { } if (this->input_type == "random_effects") { (*this->re)[i] = rnorm(this->get_expected(i), - fims_math::exp(log_sd.get_force_scalar(i))); + fims_math::exp(log_sd.get_force_scalar(i))); } if (this->input_type == "prior") { (*(this->priors[i]))[0] = diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_distribution.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_distribution.hpp index e01bb19e4..e67605f27 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_distribution.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_distribution.hpp @@ -187,7 +187,7 @@ class DnormDistributionsInterface : public DistributionsInterfaceBase { */ ParameterVector expected_mean; /** - * @brief The log_sd, which would be the standard deviation of x for the + * @brief The uncertainty, which would be the standard deviation of x for the * normal distribution. */ ParameterVector log_sd; @@ -487,7 +487,7 @@ class DnormDistributionsInterface : public DistributionsInterfaceBase { info->RegisterParameter(distribution->log_sd[i]); } if (this->log_sd[i].estimation_type_m.get() == "random_effects") { - FIMS_ERROR_LOG("log_sd values cannot be set to random effects"); + FIMS_ERROR_LOG("standard deviations cannot be set to random effects"); } } set_variable_map(&(distribution)->log_sd, this->log_sd); @@ -546,7 +546,7 @@ class DlnormDistributionsInterface : public DistributionsInterfaceBase { */ ParameterVector expected_values; /** - * @brief The log_sd, which would be the natural logarithm of the + * @brief The uncertainty, which would be the natural logarithm of the standard deviation (sd) of log(x) for this distribution. The natural log of the standard deviation is necessary because the exponential link function is applied to the log transformed standard deviation to insure diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp index f265e29bb..4ba948cc4 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp @@ -439,8 +439,7 @@ class RealVector { * @brief The constructor. */ RealVector(const RealVector& other) - : storage_m(other.storage_m), - id_m(other.id_m) {} + : storage_m(other.storage_m), id_m(other.id_m) {} /** * @brief The constructor. @@ -704,5 +703,4 @@ void set_variable_map(fims::Vector* ptr_fims_vector, std::vector> FIMSRcppInterfaceBase::fims_interface_objects; - #endif diff --git a/tests/testthat/helper-integration-tests-setup-function.R b/tests/testthat/helper-integration-tests-setup-function.R index 3b2a23dd8..158be443d 100644 --- a/tests/testthat/helper-integration-tests-setup-function.R +++ b/tests/testthat/helper-integration-tests-setup-function.R @@ -367,7 +367,6 @@ setup_and_run_FIMS_without_wrappers <- function(iter_id, recruitment_distribution$log_sd$resize(1) recruitment_distribution$log_sd[1]$value <- log(om_input[["logR_sd"]]) - recruitment_distribution$observed_values$resize(om_input[["nyr"]] - 1) recruitment_distribution$expected_values$resize(om_input[["nyr"]] - 1) for (i in 1:(om_input[["nyr"]] - 1)) { diff --git a/tests/testthat/test-rcpp-distribution.R b/tests/testthat/test-rcpp-distribution.R index 2ce85dce8..952ba46c3 100644 --- a/tests/testthat/test-rcpp-distribution.R +++ b/tests/testthat/test-rcpp-distribution.R @@ -250,7 +250,6 @@ test_that("rcpp_distribution returns correct outputs for edge cases", { dnorm_$observed_values[1]$value <- y dnorm_$expected_values[1]$value <- 0 dnorm_$log_sd[1]$value <- -10 - # evaluate the density and compare with R #' @description Test extreme log_sd values for dnorm (-10) return expected output. expect_equal(dnorm_$evaluate(), stats::dnorm(y, 0, exp(-10), TRUE)) From fcdc1e2a91b465820e54f06a5e96512f0d224d4a Mon Sep 17 00:00:00 2001 From: Andrea-Havron-NOAA <85530309+Andrea-Havron-NOAA@users.noreply.github.com> Date: Sat, 30 May 2026 00:33:47 +0000 Subject: [PATCH 4/4] chore: format feature branch --- R/fimsfit.R | 6 +- inst/include/common/fims_transformations.hpp | 549 +++++++++--------- inst/include/common/information.hpp | 84 +-- inst/include/common/model.hpp | 10 +- inst/include/common/types.hpp | 106 ++-- .../functors/density_components_base.hpp | 18 +- .../include/interface/rcpp/rcpp_interface.hpp | 213 ++++--- .../rcpp/rcpp_objects/rcpp_distribution.hpp | 8 +- .../rcpp/rcpp_objects/rcpp_interface_base.hpp | 45 +- .../rcpp/rcpp_objects/rcpp_population.hpp | 3 +- .../rcpp/rcpp_objects/rcpp_recruitment.hpp | 16 +- .../rcpp/rcpp_objects/rcpp_selectivity.hpp | 6 +- src/fims_modules.hpp | 3 +- tests/testthat/test-fimsfit.R | 4 +- 14 files changed, 526 insertions(+), 545 deletions(-) diff --git a/R/fimsfit.R b/R/fimsfit.R index f9c7e1b46..65e0195c3 100644 --- a/R/fimsfit.R +++ b/R/fimsfit.R @@ -574,7 +574,7 @@ fit_fims <- function(input, if (is.null(opt)) { failed_nlminb_object <- return_failed_nlminb(obj) - failed_fit <- fit <- FIMSFit( + failed_fit <- fit <- FIMSFit( input = input, obj = obj, opt = failed_nlminb_object[["opt"]], @@ -608,7 +608,7 @@ fit_fims <- function(input, "i" = "The failed results are being returned." )) failed_nlminb_object <- return_failed_nlminb(obj) - failed_fit <- fit <- FIMSFit( + failed_fit <- fit <- FIMSFit( input = input, obj = obj, opt = failed_nlminb_object[["opt"]], @@ -721,7 +721,7 @@ try_nlminb <- function(object, control_list, starting_values) { } return_failed_nlminb <- function(object) { - failed_nlminb_message <- c( + failed_nlminb_message <- c( "x" = "{.fun fit_fims} did not lead to a converged model.", "i" = "The resulting coefficients, probability values, or predictions are not accurate or stable and should not be used for management.", diff --git a/inst/include/common/fims_transformations.hpp b/inst/include/common/fims_transformations.hpp index 8fe276f7a..ba5655274 100644 --- a/inst/include/common/fims_transformations.hpp +++ b/inst/include/common/fims_transformations.hpp @@ -1,13 +1,13 @@ /** * @file fims_transformations.hpp * @brief Defines transformations for parameters in FIMS. - - * + + * * @details This file provides functions for applying forward and inverse * transformations to parameters, transforming parameters between input * and prior spaces, and computing log Jacobian adjustments for * change-of-variables corrections in MCMC sampling. - * + * * The transformation system supports the following transformations: * - `identity`: No transformation, parameter is on the natural scale. * - `exp`: Exponential transformation, parameter is on the log scale. @@ -15,7 +15,7 @@ * - `logit`: Logit transformation with optional bounds. * - `square`: Square transformation. * - `sqrt`: Square root transformation. - * + * * @see fims::Transformation * @see fims_math.hpp */ @@ -28,313 +28,296 @@ namespace fims_transformations { - /** - * @brief Applies a forward transformation to a scalar value. - * - * @details Maps a value from natural scale to the transformed scale - * specified by the transformation argument. The forward transformations - * and are: - * - * | Label | Forward | - * |----------|------------------| - * | identity | x | - * | log | log(x) | - * | exp | exp(x) | - * | logit | logit(x, lo, hi) | - * | square | x^2 | - * | sqrt | sqrt(x) | - * - * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). - * @param value The input value on the natural scale. - * @param transformation A fims::Transformation struct specifying the - * transformation label and any required arguments (e.g. logit bounds). - * @return The transformed value. - * @throws std::invalid_argument if the transformation label is not supported. - */ - template - inline Type ApplyTransformation( - const Type& value, - const fims::Transformation& transformation) - { - - const auto label = transformation.label; - const auto args = transformation.args; - - switch(label) { - case fims::Transformation::Label::identity: - return value; - case fims::Transformation::Label::exp: - return exp(value); - case fims::Transformation::Label::log: - return log(value); - case fims::Transformation::Label::logit: - return fims_math::logit( - Type(args.lower), - Type(args.upper), value); - case fims::Transformation::Label::square: - return value * value; - case fims::Transformation::Label::sqrt: - return sqrt(value); - // Add more cases as needed - default: - throw std::invalid_argument("Unknown transformation label"); - } +/** + * @brief Applies a forward transformation to a scalar value. + * + * @details Maps a value from natural scale to the transformed scale + * specified by the transformation argument. The forward transformations + * and are: + * + * | Label | Forward | + * |----------|------------------| + * | identity | x | + * | log | log(x) | + * | exp | exp(x) | + * | logit | logit(x, lo, hi) | + * | square | x^2 | + * | sqrt | sqrt(x) | + * + * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). + * @param value The input value on the natural scale. + * @param transformation A fims::Transformation struct specifying the + * transformation label and any required arguments (e.g. logit bounds). + * @return The transformed value. + * @throws std::invalid_argument if the transformation label is not supported. + */ +template +inline Type ApplyTransformation(const Type& value, + const fims::Transformation& transformation) { + const auto label = transformation.label; + const auto args = transformation.args; + + switch (label) { + case fims::Transformation::Label::identity: + return value; + case fims::Transformation::Label::exp: + return exp(value); + case fims::Transformation::Label::log: + return log(value); + case fims::Transformation::Label::logit: + return fims_math::logit(Type(args.lower), Type(args.upper), value); + case fims::Transformation::Label::square: + return value * value; + case fims::Transformation::Label::sqrt: + return sqrt(value); + // Add more cases as needed + default: + throw std::invalid_argument("Unknown transformation label"); } +} - - /** - * @brief Applies an inverse (back) transformation to a scalar value. - * - * @details Maps a value from the transformed scale back to the natural - * scale. This is the inverse of ApplyTransformation(). For example, if - * the input transformation is `log`, then ApplyBackTransformation applies - * `exp` to recover the natural-scale value. - * - * | Label | Back Transformation | - * |----------|----------------------| - * | identity | x | - * | log | exp(x) | - * | exp | log(x) | - * | logit | inv_logit(x, lo, hi) | - * | square | sqrt(x) | - * | sqrt | x^2 | - * - * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). - * @param value The input value on the transformed scale. - * @param transformation A fims::Transformation struct specifying the - * transformation label and any required arguments (e.g. logit bounds). - * @return The back-transformed value on the natural scale. - * @throws std::invalid_argument if the transformation label is not supported. - */ - template - inline Type ApplyBackTransformation( - const Type& value, - const fims::Transformation& transformation) { - - const auto label = transformation.label; - const auto args = transformation.args; - - Type transformed_value; - switch (label) { - case fims::Transformation::Label::identity: - transformed_value = value; - break; - case fims::Transformation::Label::exp: - transformed_value = fims_math::log(value); - break; - case fims::Transformation::Label::log: - transformed_value = fims_math::exp(value); - break; - case fims::Transformation::Label::logit: - transformed_value = fims_math::inv_logit(Type(args.lower), - Type(args.upper), value); - break; - case fims::Transformation::Label::square: - transformed_value = fims_math::sqrt(value); - break; - case fims::Transformation::Label::sqrt: - transformed_value = value * value; - break; - default: - throw std::invalid_argument( - std::string("Unknown transformation applied to a parameter, ") + - TransformationLabelToString(label) + - std::string(". Valid transformations are identity, exp, log, logit, square, and sqrt.")); - break; - } - - return transformed_value; +/** + * @brief Applies an inverse (back) transformation to a scalar value. + * + * @details Maps a value from the transformed scale back to the natural + * scale. This is the inverse of ApplyTransformation(). For example, if + * the input transformation is `log`, then ApplyBackTransformation applies + * `exp` to recover the natural-scale value. + * + * | Label | Back Transformation | + * |----------|----------------------| + * | identity | x | + * | log | exp(x) | + * | exp | log(x) | + * | logit | inv_logit(x, lo, hi) | + * | square | sqrt(x) | + * | sqrt | x^2 | + * + * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). + * @param value The input value on the transformed scale. + * @param transformation A fims::Transformation struct specifying the + * transformation label and any required arguments (e.g. logit bounds). + * @return The back-transformed value on the natural scale. + * @throws std::invalid_argument if the transformation label is not supported. + */ +template +inline Type ApplyBackTransformation( + const Type& value, const fims::Transformation& transformation) { + const auto label = transformation.label; + const auto args = transformation.args; + + Type transformed_value; + switch (label) { + case fims::Transformation::Label::identity: + transformed_value = value; + break; + case fims::Transformation::Label::exp: + transformed_value = fims_math::log(value); + break; + case fims::Transformation::Label::log: + transformed_value = fims_math::exp(value); + break; + case fims::Transformation::Label::logit: + transformed_value = + fims_math::inv_logit(Type(args.lower), Type(args.upper), value); + break; + case fims::Transformation::Label::square: + transformed_value = fims_math::sqrt(value); + break; + case fims::Transformation::Label::sqrt: + transformed_value = value * value; + break; + default: + throw std::invalid_argument( + std::string("Unknown transformation applied to a parameter, ") + + TransformationLabelToString(label) + + std::string(". Valid transformations are identity, exp, log, logit, " + "square, and sqrt.")); + break; } + return transformed_value; +} - /** - * @brief Transforms a vector of parameters from input space to prior space. - * - * @details Applies the composite transformation: - * input space -> natural scale -> prior space - * - * This is used when evaluating a prior distribution where the parameter - * is estimated in one space (e.g. log scale) but the prior is defined in - * a different space (e.g. natural scale or variance scale). - * - * If input and prior transformations are identical, the input vector is - * returned unchanged. - * - * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). - * @param input_value A vector of parameter values in input space. - * @param input The transformation applied to the parameter in input space - * (e.g. log if the parameter is stored as log(sd)). - * @param prior The transformation applied to the parameter in prior space - * (e.g. square if the prior is on variance = sd^2). - * @return A vector of parameter values transformed to prior space. - */ - template - inline fims::Vector TransformPrior(const fims::Vector& input_value, - fims::Transformation input, - fims::Transformation prior) { - - if (input.label == prior.label) { - return input_value; - } else { - - size_t n = input_value.size(); - - fims::Vector natural_parameter; - fims::Vector prior_parameter; - natural_parameter.resize(n); - prior_parameter.resize(n); +/** + * @brief Transforms a vector of parameters from input space to prior space. + * + * @details Applies the composite transformation: + * input space -> natural scale -> prior space + * + * This is used when evaluating a prior distribution where the parameter + * is estimated in one space (e.g. log scale) but the prior is defined in + * a different space (e.g. natural scale or variance scale). + * + * If input and prior transformations are identical, the input vector is + * returned unchanged. + * + * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). + * @param input_value A vector of parameter values in input space. + * @param input The transformation applied to the parameter in input space + * (e.g. log if the parameter is stored as log(sd)). + * @param prior The transformation applied to the parameter in prior space + * (e.g. square if the prior is on variance = sd^2). + * @return A vector of parameter values transformed to prior space. + */ +template +inline fims::Vector TransformPrior(const fims::Vector& input_value, + fims::Transformation input, + fims::Transformation prior) { + if (input.label == prior.label) { + return input_value; + } else { + size_t n = input_value.size(); - for(size_t i = 0; i < n; i ++) { - natural_parameter[i] = - ApplyBackTransformation(input_value[i], input); - prior_parameter[i] = - ApplyTransformation(natural_parameter[i], prior); - } + fims::Vector natural_parameter; + fims::Vector prior_parameter; + natural_parameter.resize(n); + prior_parameter.resize(n); - return prior_parameter; + for (size_t i = 0; i < n; i++) { + natural_parameter[i] = ApplyBackTransformation(input_value[i], input); + prior_parameter[i] = ApplyTransformation(natural_parameter[i], prior); } - } - - /** - * @brief Transforms a scalar parameter from input space to prior space. - * - * @details Scalar overload of TransformPrior() for use in get_observed() - * where a single element is needed. Applies the composite transformation: - * input space -> natural scale -> prior space - * - * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). - * @param input_value A scalar parameter value in input space. - * @param input The transformation applied to the parameter in input space. - * @param prior The transformation applied to the parameter in prior space. - * @return The scalar parameter value transformed to prior space. - */ - template - inline Type TransformPrior(const Type& input_value, - const fims::Transformation& input, - const fims::Transformation& prior) { - - Type natural_parameter = ApplyBackTransformation(input_value, input); - return ApplyTransformation(natural_parameter, - prior); + return prior_parameter; } +} - /** - * @brief Computes the log absolute determinant of the Jacobian for a - * change-of-variables correction in MCMC sampling. - * - * @details When a parameter is estimated in one space (input space) but - * the prior is placed in a different space (prior space), MCMC sampling - * requires a change-of-variables correction. This function computes - * log |det J| where J is the Jacobian of the composite transformation - * from input space to prior space. - * - * The Jacobian is computed using TMBad's sparse Jacobian functionality - * via SpJacFun(), which avoids allocating a dense n x n matrix. For - * element-wise transformations (the common case), the Jacobian is - * diagonal and only n non-zero entries are stored. For future multivariate - * or simplex transformations, the sparse pattern is detected automatically - * during AD taping. - * - * @note This function is only available when compiled with TMB - * (TMB_MODEL defined). For non-TMB builds it throws std::invalid_argument. - * - * @note Stan handles the Jacobian for the transformation from unconstrained - * sampling space to bounded input space automatically. This function only - * handles the additional correction from input space to prior space. - * - * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). - * @param input_value A vector of parameter values in input space. - * @param input The transformation applied to the parameter in input space - * (e.g. log if the parameter is stored as log(sd)). - * @param prior The transformation applied to the parameter in prior space - * (e.g. square if the prior is on variance = sd^2). - * @return The log absolute determinant of the Jacobian, log |det J|. - * @throws std::invalid_argument if called outside a TMB model context. - * - * @example - * // Parameter estimated as log(sd), prior on variance (sd^2) - * // Jacobian adjustment: log(2) + 2*log(sd) - * fims::Transformation log_trans, square_trans; - * log_trans.label = fims::Transformation::Label::log; - * square_trans.label = fims::Transformation::Label::square; - * Type jac = AddLogJacobian(log_sd_values, log_trans, square_trans); - */ - template - inline Type AddLogJacobian(const fims::Vector& input_value, - fims::Transformation input, - fims::Transformation prior) { - - #ifdef TMB_MODEL +/** + * @brief Transforms a scalar parameter from input space to prior space. + * + * @details Scalar overload of TransformPrior() for use in get_observed() + * where a single element is needed. Applies the composite transformation: + * input space -> natural scale -> prior space + * + * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). + * @param input_value A scalar parameter value in input space. + * @param input The transformation applied to the parameter in input space. + * @param prior The transformation applied to the parameter in prior space. + * @return The scalar parameter value transformed to prior space. + */ +template +inline Type TransformPrior(const Type& input_value, + const fims::Transformation& input, + const fims::Transformation& prior) { + Type natural_parameter = ApplyBackTransformation(input_value, input); + return ApplyTransformation(natural_parameter, prior); +} - size_t n = input_value.size(); - // Step 1: Extract double values for taping - std::vector x_val(n); - for (size_t i = 0; i < n; i++) x_val[i] = asDouble(input_value[i]); +/** + * @brief Computes the log absolute determinant of the Jacobian for a + * change-of-variables correction in MCMC sampling. + * + * @details When a parameter is estimated in one space (input space) but + * the prior is placed in a different space (prior space), MCMC sampling + * requires a change-of-variables correction. This function computes + * log |det J| where J is the Jacobian of the composite transformation + * from input space to prior space. + * + * The Jacobian is computed using TMBad's sparse Jacobian functionality + * via SpJacFun(), which avoids allocating a dense n x n matrix. For + * element-wise transformations (the common case), the Jacobian is + * diagonal and only n non-zero entries are stored. For future multivariate + * or simplex transformations, the sparse pattern is detected automatically + * during AD taping. + * + * @note This function is only available when compiled with TMB + * (TMB_MODEL defined). For non-TMB builds it throws std::invalid_argument. + * + * @note Stan handles the Jacobian for the transformation from unconstrained + * sampling space to bounded input space automatically. This function only + * handles the additional correction from input space to prior space. + * + * @tparam Type The numeric type (e.g. double, TMBad::ad_aug). + * @param input_value A vector of parameter values in input space. + * @param input The transformation applied to the parameter in input space + * (e.g. log if the parameter is stored as log(sd)). + * @param prior The transformation applied to the parameter in prior space + * (e.g. square if the prior is on variance = sd^2). + * @return The log absolute determinant of the Jacobian, log |det J|. + * @throws std::invalid_argument if called outside a TMB model context. + * + * @example + * // Parameter estimated as log(sd), prior on variance (sd^2) + * // Jacobian adjustment: log(2) + 2*log(sd) + * fims::Transformation log_trans, square_trans; + * log_trans.label = fims::Transformation::Label::log; + * square_trans.label = fims::Transformation::Label::square; + * Type jac = AddLogJacobian(log_sd_values, log_trans, square_trans); + */ +template +inline Type AddLogJacobian(const fims::Vector& input_value, + fims::Transformation input, + fims::Transformation prior) { +#ifdef TMB_MODEL + + size_t n = input_value.size(); + // Step 1: Extract double values for taping + std::vector x_val(n); + for (size_t i = 0; i < n; i++) x_val[i] = asDouble(input_value[i]); + + // Step 2: Tape the composite transformation into an ADFun + TMBad::ADFun<> f; + { + // Start recording - // Step 2: Tape the composite transformation into an ADFun - TMBad::ADFun<> f; - { - // Start recording - std::vector x_ad(n); for (size_t i = 0; i < n; i++) x_ad[i] = x_val[i]; - + f.glob.ad_start(); f.glob.Independent(x_ad); // vector version, not scalar - + std::vector y_ad(n); for (size_t i = 0; i < n; i++) { - TMBad::ad_aug natural = ApplyBackTransformation(x_ad[i], input); - y_ad[i] = ApplyTransformation(natural, prior); + TMBad::ad_aug natural = ApplyBackTransformation(x_ad[i], input); + y_ad[i] = ApplyTransformation(natural, prior); } - + TMBad::Dependent(y_ad); // vector version f.glob.ad_stop(); - } - - // Step 3: Generate sparse Jacobian function - TMBad::Sparse> spjac = f.SpJacFun(); + } - // Step 4: Evaluate sparse Jacobian at current values - std::vector jac_vals = spjac(x_val); + // Step 3: Generate sparse Jacobian function + TMBad::Sparse> spjac = f.SpJacFun(); - // Step 5: Compute log abs det from sparse values - // For diagonal case, all non-zeros are diagonal entries - // For general case, need to use spjac.i and spjac.j for row/col indices - Type log_abs_det = Type(0.0); + // Step 4: Evaluate sparse Jacobian at current values + std::vector jac_vals = spjac(x_val); - if (n == spjac.i.size()) { - // Diagonal case - all non-zeros are on diagonal - for (size_t k = 0; k < jac_vals.size(); k++) { - log_abs_det += fims_math::log(fims_math::ad_fabs(Type(jac_vals[k]))); - } - } else { - // General case - need log determinant of sparse matrix - // Build Eigen sparse matrix from spjac pattern - Eigen::SparseMatrix J(n, n); - std::vector> triplets; - for (size_t k = 0; k < jac_vals.size(); k++) { - triplets.push_back(Eigen::Triplet( - spjac.i[k], spjac.j[k], jac_vals[k])); - } - J.setFromTriplets(triplets.begin(), triplets.end()); + // Step 5: Compute log abs det from sparse values + // For diagonal case, all non-zeros are diagonal entries + // For general case, need to use spjac.i and spjac.j for row/col indices + Type log_abs_det = Type(0.0); - // For sparse log determinant use Eigen's sparse LU - Eigen::SparseLU> solver; - solver.compute(J); - log_abs_det = Type(solver.logAbsDeterminant()); + if (n == spjac.i.size()) { + // Diagonal case - all non-zeros are on diagonal + for (size_t k = 0; k < jac_vals.size(); k++) { + log_abs_det += fims_math::log(fims_math::ad_fabs(Type(jac_vals[k]))); } + } else { + // General case - need log determinant of sparse matrix + // Build Eigen sparse matrix from spjac pattern + Eigen::SparseMatrix J(n, n); + std::vector> triplets; + for (size_t k = 0; k < jac_vals.size(); k++) { + triplets.push_back( + Eigen::Triplet(spjac.i[k], spjac.j[k], jac_vals[k])); + } + J.setFromTriplets(triplets.begin(), triplets.end()); - return log_abs_det; - - #else - - throw std::invalid_argument( - std::string("Jacobian adjustments currently only work for TMB models.") - ); - #endif //TMB_MODEL + // For sparse log determinant use Eigen's sparse LU + Eigen::SparseLU> solver; + solver.compute(J); + log_abs_det = Type(solver.logAbsDeterminant()); } -} // namespace fims_transformations + return log_abs_det; + +#else + + throw std::invalid_argument( + std::string("Jacobian adjustments currently only work for TMB models.")); +#endif // TMB_MODEL +} + +} // namespace fims_transformations #endif /* FIMS_COMMON_FIMS_TRANSFORMATIONS_HPP */ \ No newline at end of file diff --git a/inst/include/common/information.hpp b/inst/include/common/information.hpp index fd6f43107..153450b68 100644 --- a/inst/include/common/information.hpp +++ b/inst/include/common/information.hpp @@ -138,50 +138,50 @@ class Information { typedef typename std::unordered_map< uint32_t, std::shared_ptr>>::iterator model_map_iterator; /**< iterator for variable map>*/ - -/** - * @brief A structure to hold a pointer to a parameter vector and its transformation - * metadata for use in the variable map. - * - * @details Each entry in the variable map corresponds to a single parameter - * vector and stores a pointer to the parameter values along with two - * transformation labels: - * - `input_transformation`: the transformation applied to the parameter - * in the input space (e.g. log, logit). This is the space in which - * the parameter is estimated. - * - `prior_transformation`: the transformation applied to the parameter - * in the prior space (e.g. identity, square). This is the space in - * which the prior distribution is defined. - */ -struct VariableMapEntry { - fims::Vector* variable = nullptr; - fims::Transformation input_transformation; - fims::Transformation prior_transformation; -/** - * @brief Constructor for VariableMapEntry. - */ -VariableMapEntry() { - input_transformation.label = fims::Transformation::Label::log; - prior_transformation.label = fims::Transformation::Label::log; -} + /** + * @brief A structure to hold a pointer to a parameter vector and its + * transformation metadata for use in the variable map. + * + * @details Each entry in the variable map corresponds to a single parameter + * vector and stores a pointer to the parameter values along with two + * transformation labels: + * - `input_transformation`: the transformation applied to the parameter + * in the input space (e.g. log, logit). This is the space in which + * the parameter is estimated. + * - `prior_transformation`: the transformation applied to the parameter + * in the prior space (e.g. identity, square). This is the space in + * which the prior distribution is defined. + */ + struct VariableMapEntry { + fims::Vector* variable = nullptr; + fims::Transformation input_transformation; + fims::Transformation prior_transformation; + + /** + * @brief Constructor for VariableMapEntry. + */ + VariableMapEntry() { + input_transformation.label = fims::Transformation::Label::log; + prior_transformation.label = fims::Transformation::Label::log; + } -/** - * @brief Constructor for VariableMapEntry with all fields initialized. - * - * @param variable Pointer to the fims::Vector holding the parameter values. - * @param input_transformation The transformation applied to the parameter - * in the input space (e.g. log, logit). - * @param prior_transformation The transformation applied to the parameter - * in the prior space (e.g. identity, square). - */ - VariableMapEntry(fims::Vector* variable, - fims::Transformation input_transformation, - fims::Transformation prior_transformation) - : variable(variable), - input_transformation(input_transformation), - prior_transformation(prior_transformation) {} -}; + /** + * @brief Constructor for VariableMapEntry with all fields initialized. + * + * @param variable Pointer to the fims::Vector holding the parameter values. + * @param input_transformation The transformation applied to the parameter + * in the input space (e.g. log, logit). + * @param prior_transformation The transformation applied to the parameter + * in the prior space (e.g. identity, square). + */ + VariableMapEntry(fims::Vector* variable, + fims::Transformation input_transformation, + fims::Transformation prior_transformation) + : variable(variable), + input_transformation(input_transformation), + prior_transformation(prior_transformation) {} + }; std::unordered_map variable_map; /**input_type == "prior") { nll_vec[nll_vec_idx] = -d->evaluate(); if (this->jacobian_flag) { - for(size_t i=0; i < d->key.size(); i++) { - if ((*(d->input_transformation[i])).label != + for (size_t i = 0; i < d->key.size(); i++) { + if ((*(d->input_transformation[i])).label != (*(d->prior_transformation[i])).label) { - nll_vec[nll_vec_idx] += - -fims_transformations::AddLogJacobian( - *(d->priors[i]), - *(d->input_transformation[i]), + nll_vec[nll_vec_idx] += -fims_transformations::AddLogJacobian( + *(d->priors[i]), *(d->input_transformation[i]), *(d->prior_transformation[i])); } } diff --git a/inst/include/common/types.hpp b/inst/include/common/types.hpp index 737851db5..c1ad70f44 100644 --- a/inst/include/common/types.hpp +++ b/inst/include/common/types.hpp @@ -1,8 +1,9 @@ /** * @file types.hpp - * @brief Centralized definition of core model types, enumerations, and structures. - * @details This file serves as a lightweight, dependency-free foundation for the - * fims namespace. + * @brief Centralized definition of core model types, enumerations, and + * structures. + * @details This file serves as a lightweight, dependency-free foundation for + * the fims namespace. * @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. @@ -15,9 +16,9 @@ namespace fims { /** * @brief A data structure for transforming parameters - * - * @details The `Transformation` struct defines the type of transformation to - * apply to a parameter and any necessary arguments for that transformation. + * + * @details The `Transformation` struct defines the type of transformation to + * apply to a parameter and any necessary arguments for that transformation. * This allows for flexible parameter transformations within the FIMS framework. */ struct Transformation { @@ -25,51 +26,51 @@ struct Transformation { * @brief An enumeration of supported transformation types. * @details Supported transformations include: * - `identity`: No transformation, parameter is on the natural scale. - * - `exp`: Exponential transformation, parameter is on the multiplicative scale. + * - `exp`: Exponential transformation, parameter is on the multiplicative + * scale. * - `log`: Log transformation, parameter is on the log scale. * - `logit`: Logit transformation, parameter is on the logit scale. * - `square`: Square transformation, parameter is on the squared scale. - * - `sqrt`: Square root transformation, parameter is on the square root scale. + * - `sqrt`: Square root transformation, parameter is on the square root + * scale. */ - enum class Label { - identity, - exp, - log, - logit, - square, - sqrt - }; - /** - * @brief A structure to hold arguments for transformations that require them. - * @details For the logit transformation, `lower` and `upper` specify the - * bounds of the parameter in natural space. - */ - struct Args { - //used by logit - double lower = 0.0; - double upper = 1.0; - }; - - - Label label = Label::log; - Args args{}; - + enum class Label { identity, exp, log, logit, square, sqrt }; + /** + * @brief A structure to hold arguments for transformations that require them. + * @details For the logit transformation, `lower` and `upper` specify the + * bounds of the parameter in natural space. + */ + struct Args { + // used by logit + double lower = 0.0; + double upper = 1.0; }; - static const char* TransformationLabelToString( - fims::Transformation::Label label) { - switch (label) { - case fims::Transformation::Label::identity: return "identity"; - case fims::Transformation::Label::exp: return "exp"; - case fims::Transformation::Label::log: return "log"; - case fims::Transformation::Label::logit: return "logit"; - case fims::Transformation::Label::square: return "square"; - case fims::Transformation::Label::sqrt: return "sqrt"; - default: return "unknown"; - } + Label label = Label::log; + Args args{}; +}; + +static const char* TransformationLabelToString( + fims::Transformation::Label label) { + switch (label) { + case fims::Transformation::Label::identity: + return "identity"; + case fims::Transformation::Label::exp: + return "exp"; + case fims::Transformation::Label::log: + return "log"; + case fims::Transformation::Label::logit: + return "logit"; + case fims::Transformation::Label::square: + return "square"; + case fims::Transformation::Label::sqrt: + return "sqrt"; + default: + return "unknown"; } +} - /** +/** * @brief A structure to hold distribution type information for priors. * @details Supported distributions include: * - `Normal`: Normal distribution. @@ -82,24 +83,19 @@ struct Distribution { /** * @brief An enumeration of supported distribution types. */ - enum class Label { - Normal, - Lognormal, - Gamma, - InvGamma - }; + enum class Label { Normal, Lognormal, Gamma, InvGamma }; Label label = Label::Normal; }; fims::Distribution::Label StringToDistributionLabel(const std::string& name) { - if (name == "Normal") return fims::Distribution::Label::Normal; - if (name == "Lognormal") return fims::Distribution::Label::Lognormal; - if (name == "Gamma") return fims::Distribution::Label::Gamma; - if (name == "InvGamma") return fims::Distribution::Label::InvGamma; - throw std::invalid_argument("Unsupported distribution: " + name); + if (name == "Normal") return fims::Distribution::Label::Normal; + if (name == "Lognormal") return fims::Distribution::Label::Lognormal; + if (name == "Gamma") return fims::Distribution::Label::Gamma; + if (name == "InvGamma") return fims::Distribution::Label::InvGamma; + throw std::invalid_argument("Unsupported distribution: " + name); } -} // namespace fims +} // namespace fims #endif /* FIMS_COMMON_TYPES_HPP */ \ No newline at end of file diff --git a/inst/include/distributions/functors/density_components_base.hpp b/inst/include/distributions/functors/density_components_base.hpp index 5bda02799..fa364f14f 100644 --- a/inst/include/distributions/functors/density_components_base.hpp +++ b/inst/include/distributions/functors/density_components_base.hpp @@ -59,7 +59,7 @@ struct DensityComponentBase : public fims_model_object::FIMSObject { /** @brief Pointer to input transformation */ std::vector input_transformation; - + /** @brief Pointer to prior transformation */ std::vector prior_transformation; @@ -88,22 +88,20 @@ struct DensityComponentBase : public fims_model_object::FIMSObject { throw std::runtime_error("No priors defined for this distribution."); } else if (priors.size() == 1) { this->transformed_priors.resize((*priors[0]).size()); - this->transformed_priors = - fims_transformations::TransformPrior(*(priors[0]), - *(input_transformation[0]), - *(prior_transformation[0])); + this->transformed_priors = fims_transformations::TransformPrior( + *(priors[0]), *(input_transformation[0]), + *(prior_transformation[0])); } else if (priors.size() > 1) { size_t n = priors.size(); this->transformed_priors.resize(n); for (size_t i = 0; i < n; i++) { - transformed_priors[i] = - fims_transformations::TransformPrior((*(priors[i]))[0], - *(input_transformation[i]), - *(prior_transformation[i])); + transformed_priors[i] = fims_transformations::TransformPrior( + (*(priors[i]))[0], *(input_transformation[i]), + *(prior_transformation[i])); + } } } } -} /** * @brief Retrieve one observed value based on `input_type`. diff --git a/inst/include/interface/rcpp/rcpp_interface.hpp b/inst/include/interface/rcpp/rcpp_interface.hpp index 4863703b2..6417421ad 100644 --- a/inst/include/interface/rcpp/rcpp_interface.hpp +++ b/inst/include/interface/rcpp/rcpp_interface.hpp @@ -96,11 +96,11 @@ bool CreateTMBModel() { info0->CheckModel(); info->CreateModel(); -/* - // instantiate the model? TODO: Ask Matthew what this does - std::shared_ptr> m0 = - fims_model::Model::GetInstance(); -*/ + /* + // instantiate the model? TODO: Ask Matthew what this does + std::shared_ptr> m0 = + fims_model::Model::GetInstance(); + */ return true; } @@ -229,11 +229,13 @@ Rcpp::List get_random_names(Rcpp::List pars) { } /** - * @brief Classifies an R formula's Left-Hand Side (LHS) and extracts the transformation. - * * @details This helper function inspects the Abstract Syntax Tree (AST) of the raw - * R expression pointer (SEXP) passed from the formula's left-hand side. It maps - * recognized R mathematical operations (like `log()`, `exp()`, or exponentiation `^2`) - * to their respective domain-specific strongly typed enums defined in the fims namespace. + * @brief Classifies an R formula's Left-Hand Side (LHS) and extracts the + * transformation. + * * @details This helper function inspects the Abstract Syntax Tree (AST) of + * the raw R expression pointer (SEXP) passed from the formula's left-hand side. + * It maps recognized R mathematical operations (like `log()`, `exp()`, or + * exponentiation `^2`) to their respective domain-specific strongly typed enums + * defined in the fims namespace. * * Supported AST structural mappings include: * - \b Symbol (e.g., \code{sd}): Maps to \c Label::identity * - \b Language Call (e.g., \code{log(sd)}): Maps to \c Label::log @@ -242,107 +244,104 @@ Rcpp::List get_random_names(Rcpp::List pars) { * - \b Language Call (e.g., \code{sqrt(sd)}): Maps to \c Label::sqrt * - \b Infix Exponent (e.g., \code{sd^2}): Maps to \c Label::square * - * @param lhs_raw A raw \code{SEXP} representing the left-hand side node of an R formula. - * * @throws std::runtime_error (via \code{Rcpp::stop}) if the expression tree structure - * is unrecognizable or contains an unmapped transformation function. + * @param lhs_raw A raw \code{SEXP} representing the left-hand side node of an R + * formula. + * * @throws std::runtime_error (via \code{Rcpp::stop}) if the expression tree + * structure is unrecognizable or contains an unmapped transformation function. * * @return A \code{fims::Transformation::Label} with the enumeration token. */ fims::Transformation classify_and_extract_transformation(SEXP lhs_raw) { - fims::Transformation trans_config; - - if (Rcpp::is(lhs_raw)) { - trans_config.label = fims::Transformation::Label::identity; + + if (Rcpp::is(lhs_raw)) { + trans_config.label = fims::Transformation::Label::identity; + return trans_config; + } + + if (Rcpp::is(lhs_raw)) { + Rcpp::Language lhs_lang = Rcpp::as(lhs_raw); + Rcpp::CharacterVector op_name = + Rcpp::as(Rcpp::as(lhs_lang[0])); + + // Match R operators/functions to your fims::Transformation::Label enums + if (op_name[0] == "log") { + trans_config.label = fims::Transformation::Label::log; + return trans_config; + } else if (op_name[0] == "exp") { + trans_config.label = fims::Transformation::Label::exp; + return trans_config; + } else if (op_name[0] == "logit") { + trans_config.label = fims::Transformation::Label::logit; + // Extract bounds if provided, otherwise default to 0 and 1 + if (lhs_lang.size() >= 3) { + trans_config.args.lower = Rcpp::as(lhs_lang[2]); + } + if (lhs_lang.size() >= 4) { + trans_config.args.upper = Rcpp::as(lhs_lang[3]); + } + // Validate bounds + if (trans_config.args.lower >= trans_config.args.upper) { + Rcpp::stop("Logit transformation requires lower < upper."); + } + return trans_config; + } else if (op_name[0] == "sqrt") { + trans_config.label = fims::Transformation::Label::sqrt; + return trans_config; + } else if (op_name[0] == "^") { + int power = Rcpp::as(lhs_lang[2]); + if (power == 2) { + trans_config.label = fims::Transformation::Label::square; return trans_config; + } } - - if (Rcpp::is(lhs_raw)) { - Rcpp::Language lhs_lang = Rcpp::as(lhs_raw); - Rcpp::CharacterVector op_name = Rcpp::as(Rcpp::as(lhs_lang[0])); - - // Match R operators/functions to your fims::Transformation::Label enums - if (op_name[0] == "log") { - trans_config.label = fims::Transformation::Label::log; - return trans_config; - } - else if (op_name[0] == "exp") { - trans_config.label = fims::Transformation::Label::exp; - return trans_config; - } - else if (op_name[0] == "logit") { - trans_config.label = fims::Transformation::Label::logit; - // Extract bounds if provided, otherwise default to 0 and 1 - if (lhs_lang.size() >= 3) { - trans_config.args.lower = Rcpp::as(lhs_lang[2]); - } - if (lhs_lang.size() >= 4) { - trans_config.args.upper = Rcpp::as(lhs_lang[3]); - } - // Validate bounds - if (trans_config.args.lower >= trans_config.args.upper) { - Rcpp::stop("Logit transformation requires lower < upper."); - } - return trans_config; - } - else if (op_name[0] == "sqrt") { - trans_config.label = fims::Transformation::Label::sqrt; - return trans_config; - } - else if (op_name[0] == "^") { - int power = Rcpp::as(lhs_lang[2]); - if (power == 2) { - trans_config.label = fims::Transformation::Label::square; - return trans_config; - } - } - // Add reciprocal parsing here if you still want to route 1/x variants - } - - // Fallback/Error state if formula layout isn't supported - Rcpp::stop("Unsupported Left-Hand Side (LHS) transformation operator."); -} + // Add reciprocal parsing here if you still want to route 1/x variants + } + // Fallback/Error state if formula layout isn't supported + Rcpp::stop("Unsupported Left-Hand Side (LHS) transformation operator."); +} /** * @brief A structure to hold the parsed components of a distribution formula. */ struct FormulaComponents { - fims::Transformation transformation; - fims::Distribution::Label distribution; - std::vector hyperparameters; + fims::Transformation transformation; + fims::Distribution::Label distribution; + std::vector hyperparameters; }; //' Parse a Distributional Formula //' -//' Parses an R formula specifying a target variable and its prior/likelihood distribution -//' along with its parameters (e.g., \code{y ~ dnorm(0, 1)}). +//' Parses an R formula specifying a target variable and its prior/likelihood +//distribution ' along with its parameters (e.g., \code{y ~ dnorm(0, 1)}). //' -//' @param f A standard R \code{Formula} object. It must follow the structure +//' @param f A standard R \code{Formula} object. It must follow the structure //' \code{variable ~ distribution(param1, param2, ...)}. //' //' @details -//' The function unpacks the R formula by treating it as an abstract syntax tree (AST) -//' via the \code{Rcpp::Language} class: -//' \itemize{ -//' \item \strong{Operator (\code{[0]}):} The tilde (\code{~}) operator. -//' \item \strong{LHS (\code{[1]}):} Extracted as a \code{Symbol} and converted to a character vector representing the target variable name. -//' \item \strong{RHS (\code{[2]}):} Treated as a nested \code{Language} call where the head (\code{[0]}) is the distribution name, and subsequent elements are the numeric parameters. -//' } +//' The function unpacks the R formula by treating it as an abstract syntax tree +//(AST) ' via the \code{Rcpp::Language} class: ' \itemize{ ' \item +//\strong{Operator (\code{[0]}):} The tilde (\code{~}) operator. ' \item +//\strong{LHS (\code{[1]}):} Extracted as a \code{Symbol} and converted to a +//character vector representing the target variable name. ' \item \strong{RHS +//(\code{[2]}):} Treated as a nested \code{Language} call where the head +//(\code{[0]}) is the distribution name, and subsequent elements are the numeric +//parameters. ' } //' //' @return A named \code{Rcpp::List} containing three elements: //' \itemize{ -//' \item \code{variable}: A character vector holding the name of the LHS variable. -//' \item \code{distribution}: A character vector holding the name of the RHS distribution function. -//' \item \code{hyperparameters}: A numeric vector containing the extracted hyperparameter values. -//' } +//' \item \code{variable}: A character vector holding the name of the LHS +//variable. ' \item \code{distribution}: A character vector holding the name +//of the RHS distribution function. ' \item \code{hyperparameters}: A numeric +//vector containing the extracted hyperparameter values. ' } FormulaComponents parse_distribution_formula(Rcpp::Formula f) { - // Convert the Formula to a standard standard R language object (Call) Rcpp::Language formula = Rcpp::as(f); // Parse the formula // R formulas are structured as a tree where element 0 is the operator `~` - // element 1 is the Left-Hand Side (LHS), and element 2 is the Right-Hand Side (RHS) + // element 1 is the Left-Hand Side (LHS), and element 2 is the Right-Hand Side + // (RHS) SEXP lhs_raw = formula[1]; @@ -351,40 +350,40 @@ FormulaComponents parse_distribution_formula(Rcpp::Formula f) { // 2. Parse Right-Hand Side (RHS) if (!Rcpp::is(formula[2])) { - Rcpp::stop("Right-Hand Side (RHS) must be a distribution function call."); + Rcpp::stop("Right-Hand Side (RHS) must be a distribution function call."); } Rcpp::Language rhs = Rcpp::as(formula[2]); - - // The first element of a function call language object is the function name itself - // Convert formula->Symbol->CharacterVector->enum + + // The first element of a function call language object is the function name + // itself Convert formula->Symbol->CharacterVector->enum Rcpp::Symbol distribution_symbol = Rcpp::as(rhs[0]); - Rcpp::CharacterVector distribution_string = Rcpp::as(distribution_symbol); + Rcpp::CharacterVector distribution_string = + Rcpp::as(distribution_symbol); // 3. Extract the numeric parameters from the function arguments - int num_args = rhs.size() - 1; // Subtract 1 because element 0 is the function name + int num_args = + rhs.size() - 1; // Subtract 1 because element 0 is the function name std::vector hyperparameters(num_args); - - for(int i = 0; i < num_args; ++i) { + + for (int i = 0; i < num_args; ++i) { SEXP arg = rhs[i + 1]; if (Rcpp::is(arg) || Rcpp::is(arg)) { - // Evaluate the expression to get a numeric value - Rcpp::Environment base_env = Rcpp::Environment::base_env(); - SEXP result = Rcpp::Rcpp_eval(arg, base_env); - hyperparameters[i] = Rcpp::as(result); + // Evaluate the expression to get a numeric value + Rcpp::Environment base_env = Rcpp::Environment::base_env(); + SEXP result = Rcpp::Rcpp_eval(arg, base_env); + hyperparameters[i] = Rcpp::as(result); } else { - hyperparameters[i] = Rcpp::as(arg); + hyperparameters[i] = Rcpp::as(arg); } } - return FormulaComponents{ - transform, - fims::StringToDistributionLabel(Rcpp::as(distribution_string)), - hyperparameters - }; + return FormulaComponents{transform, + fims::StringToDistributionLabel( + Rcpp::as(distribution_string)), + hyperparameters}; } void setup_prior(Rcpp::Formula f, Rcpp::List parameter_vectors) { - Rcpp::IntegerVector ids(parameter_vectors.size()); // Parse formula into components @@ -397,14 +396,13 @@ void setup_prior(Rcpp::Formula f, Rcpp::List parameter_vectors) { // Extract the raw pointer from the XPtr for (int i = 0; i < parameter_vectors.size(); i++) { ParameterVector pv = Rcpp::as(parameter_vectors[i]); - + // Modifies shared Transformation object *pv.prior_transformation_m = prior_transformation; ids[i] = pv.id_m; } - - switch (distribution_name) { + switch (distribution_name) { case fims::Distribution::Label::Normal: { auto prior = std::make_shared(); prior->expected_values[0].initial_value_m = hyperparameters[0]; @@ -422,13 +420,14 @@ void setup_prior(Rcpp::Formula f, Rcpp::List parameter_vectors) { } default: - throw std::invalid_argument("Unsupported distribution type in add_prior."); + throw std::invalid_argument( + "Unsupported distribution type in add_prior."); } } -//Define setup_prior_function in ParameterVector class +// Define setup_prior_function in ParameterVector class inline void register_prior_functions() { - ParameterVector::setup_prior_function = setup_prior; + ParameterVector::setup_prior_function = setup_prior; } void add_shared_prior(Rcpp::Formula f, Rcpp::List parameter_vectors) { diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_distribution.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_distribution.hpp index e67605f27..76fab81fd 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_distribution.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_distribution.hpp @@ -507,8 +507,10 @@ class DnormDistributionsInterface : public DistributionsInterfaceBase { FIMS_ERROR_LOG("expected_mean cannot be set to random effects"); } } - expected_mean.input_transformation_m->label = fims::Transformation::Label::identity; - expected_mean.prior_transformation_m->label = fims::Transformation::Label::identity; + expected_mean.input_transformation_m->label = + fims::Transformation::Label::identity; + expected_mean.prior_transformation_m->label = + fims::Transformation::Label::identity; set_variable_map(&(distribution)->expected_mean, this->expected_mean); info->density_components[distribution->id] = distribution; @@ -830,7 +832,7 @@ class DlnormDistributionsInterface : public DistributionsInterfaceBase { FIMS_ERROR_LOG("standard deviations cannot be set to random effects"); } } - set_variable_map(&(distribution)->log_sd, this->log_sd); + set_variable_map(&(distribution)->log_sd, this->log_sd); info->density_components[distribution->id] = distribution; diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp index 4ba948cc4..fce905ce3 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp @@ -154,13 +154,13 @@ class ParameterVector { /** * @brief The transformation of the parameter in the prior distribution. */ - std::shared_ptr input_transformation_m = - std::make_shared(); + std::shared_ptr input_transformation_m = + std::make_shared(); /** * @brief The transformation of the parameter in the prior distribution. */ - std::shared_ptr prior_transformation_m = - std::make_shared(); + std::shared_ptr prior_transformation_m = + std::make_shared(); /** * @brief The local ID of the Parameter object. */ @@ -272,15 +272,15 @@ class ParameterVector { return Rcpp::wrap(this->storage_m->at(pos - 1)); } - // register function first to avoid issues with circular dependancies - // in the include chain. The setup_prior_function will be set to equal the + // register function first to avoid issues with circular dependancies + // in the include chain. The setup_prior_function will be set to equal the // setup_prior() function in rcpp_interface.hpp. static std::function setup_prior_function; void add_prior(Rcpp::Formula f) { - if (setup_prior_function) { - Rcpp::List pv_list = Rcpp::List::create(Rcpp::wrap(*this)); - setup_prior_function(f, pv_list); + if (setup_prior_function) { + Rcpp::List pv_list = Rcpp::List::create(Rcpp::wrap(*this)); + setup_prior_function(f, pv_list); } else { throw std::runtime_error("setup_prior_function not registered."); } @@ -383,9 +383,9 @@ class ParameterVector { } }; // static member definitions -uint32_t ParameterVector::id_g = 0; -std::function ParameterVector::setup_prior_function = nullptr; - +uint32_t ParameterVector::id_g = 0; +std::function + ParameterVector::setup_prior_function = nullptr; /** * @brief Output for std::ostream& for a ParameterVector. @@ -685,20 +685,17 @@ class FIMSRcppInterfaceBase { return ss.str(); } -template -void set_variable_map(fims::Vector* ptr_fims_vector, - ParameterVector parameter_vector) { + template + void set_variable_map(fims::Vector* ptr_fims_vector, + ParameterVector parameter_vector) { + std::shared_ptr> info = + fims_info::Information::GetInstance(); - std::shared_ptr> info = - fims_info::Information::GetInstance(); - - info->variable_map[parameter_vector.id_m] = - typename fims_info::Information::VariableMapEntry( - ptr_fims_vector, - *parameter_vector.input_transformation_m, - *parameter_vector.prior_transformation_m); + info->variable_map[parameter_vector.id_m] = + typename fims_info::Information::VariableMapEntry( + ptr_fims_vector, *parameter_vector.input_transformation_m, + *parameter_vector.prior_transformation_m); } - }; std::vector> FIMSRcppInterfaceBase::fims_interface_objects; diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp index 7ed5806dc..44d01ba09 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp @@ -417,7 +417,8 @@ class PopulationInterface : public PopulationInterfaceBase { info->RegisterRandomEffect(population->log_init_naa[i]); } } - info->variable_map[this->log_init_naa.id_m].variable = &(population)->log_init_naa; + info->variable_map[this->log_init_naa.id_m].variable = + &(population)->log_init_naa; for (size_t i = 0; i < ages.size(); i++) { population->ages[i] = this->ages[i]; diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp index 9251b533c..b0319459b 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp @@ -353,12 +353,14 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { info->RegisterRandomEffect(recruitment->logit_steep[i]); } } - - //set transformations for parameter since default is log - logit_steep.input_transformation_m->label = fims::Transformation::Label::logit; + + // set transformations for parameter since default is log + logit_steep.input_transformation_m->label = + fims::Transformation::Label::logit; logit_steep.input_transformation_m->args.lower = 0.2; logit_steep.input_transformation_m->args.upper = 1.0; - logit_steep.prior_transformation_m->label = fims::Transformation::Label::logit; + logit_steep.prior_transformation_m->label = + fims::Transformation::Label::logit; logit_steep.prior_transformation_m->args.lower = 0.2; logit_steep.prior_transformation_m->args.upper = 1.0; @@ -382,7 +384,8 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { info->RegisterRandomEffect(recruitment->log_rzero[i]); } } - info->variable_map[this->log_rzero.id_m].variable = &(recruitment)->log_rzero; + info->variable_map[this->log_rzero.id_m].variable = + &(recruitment)->log_rzero; // set log_recruit_devs recruitment->log_recruit_devs.resize(this->log_devs.size()); for (size_t i = 0; i < this->log_devs.size(); i++) { @@ -404,7 +407,8 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { } } - info->variable_map[this->log_devs.id_m].variable = &(recruitment)->log_recruit_devs; + info->variable_map[this->log_devs.id_m].variable = + &(recruitment)->log_recruit_devs; // set log_r recruitment->log_r.resize(this->log_r.size()); diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp index ab9c675b7..cd9ab384b 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp @@ -536,7 +536,8 @@ class DoubleLogisticSelectivityInterface : public SelectivityInterfaceBase { info->RegisterRandomEffect(selectivity->slope_asc[i]); } } - info->variable_map[this->slope_asc.id_m].variable = &(selectivity)->slope_asc; + info->variable_map[this->slope_asc.id_m].variable = + &(selectivity)->slope_asc; selectivity->inflection_point_desc.resize( this->inflection_point_desc.size()); @@ -584,7 +585,8 @@ class DoubleLogisticSelectivityInterface : public SelectivityInterfaceBase { } } - info->variable_map[this->slope_desc.id_m].variable = &(selectivity)->slope_desc; + info->variable_map[this->slope_desc.id_m].variable = + &(selectivity)->slope_desc; // add to Information info->selectivity_models[selectivity->id] = selectivity; diff --git a/src/fims_modules.hpp b/src/fims_modules.hpp index 6e25c9225..07ffb2685 100644 --- a/src/fims_modules.hpp +++ b/src/fims_modules.hpp @@ -89,7 +89,8 @@ RCPP_MODULE(fims) { "get_random_names", &get_random_names, "See " "https://noaa-fims.github.io/FIMS/doxygen/rcpp__interface_8hpp.html."); - Rcpp::function("add_shared_prior", &setup_prior, + Rcpp::function( + "add_shared_prior", &setup_prior, "See " "https://noaa-fims.github.io/FIMS/doxygen/rcpp__interface_8hpp.html."); Rcpp::function( diff --git a/tests/testthat/test-fimsfit.R b/tests/testthat/test-fimsfit.R index b27410c23..8f2439f0d 100644 --- a/tests/testthat/test-fimsfit.R +++ b/tests/testthat/test-fimsfit.R @@ -238,7 +238,7 @@ test_that("fit_fims() errors when optimization fails to converge", { value = -Inf ), by = c("module_name", "label", "age") - ) |> + ) |> initialize_fims(data = data_4_model) test_results <- suppressWarnings(suppressMessages( fit_fims(initialized_poor_model, optimize = TRUE) @@ -249,6 +249,6 @@ test_that("fit_fims() errors when optimization fails to converge", { names(get_opt(test_results)), c("par", "objective", "convergence", "message") ) - + clear() })