diff --git a/CMakeLists.txt b/CMakeLists.txt index 36f8a81ed..c6d1cbe1d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -126,7 +126,7 @@ target_include_directories(fims_test INTERFACE inst/include ${rcpp_SOURCE_DIR}/inst/include - ${R_HOME}/include + ${R_HOME}/include ) # 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/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 new file mode 100644 index 000000000..ba5655274 --- /dev/null +++ b/inst/include/common/fims_transformations.hpp @@ -0,0 +1,323 @@ +/** + * @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/information.hpp b/inst/include/common/information.hpp index 55bfb2c29..153450b68 100644 --- a/inst/include/common/information.hpp +++ b/inst/include/common/information.hpp @@ -139,10 +139,54 @@ class Information { uint32_t, std::shared_ptr>>::iterator model_map_iterator; /**< iterator for variable map>*/ - std::unordered_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 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..13ce23ad7 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,19 @@ 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..c1ad70f44 --- /dev/null +++ b/inst/include/common/types.hpp @@ -0,0 +1,101 @@ +/** + * @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..fa364f14f 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,26 @@ 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 +117,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/interface/rcpp/rcpp_interface.hpp b/inst/include/interface/rcpp/rcpp_interface.hpp index e35369c64..6417421ad 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(); - + /* + // instantiate the model? TODO: Ask Matthew what this does + std::shared_ptr> m0 = + fims_model::Model::GetInstance(); + */ return true; } @@ -226,6 +228,212 @@ 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; + } + + 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..76fab81fd 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_distribution.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_distribution.hpp @@ -490,7 +490,7 @@ class DnormDistributionsInterface : 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); distribution->use_mean = this->use_mean_m.get(); distribution->expected_mean.resize(this->expected_mean.size()); @@ -507,8 +507,11 @@ 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; @@ -829,7 +832,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..fce905ce3 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 { } } }; +// static member definitions uint32_t ParameterVector::id_g = 0; +std::function + ParameterVector::setup_prior_function = nullptr; /** * @brief Output for std::ostream& for a ParameterVector. @@ -645,6 +684,18 @@ 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; 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..44d01ba09 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,8 @@ 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..b0319459b 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp @@ -353,7 +353,16 @@ 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 +384,8 @@ 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 +407,8 @@ 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 +429,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..cd9ab384b 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,8 @@ 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 +562,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 +585,8 @@ 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..07ffb2685 100644 --- a/src/fims_modules.hpp +++ b/src/fims_modules.hpp @@ -89,6 +89,10 @@ 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 +169,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/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() }) 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()