feat(edm): implement EDM prediction framework (Simplex, S-Map, GP-EDM)

.Rbuildignore

@@ -22,3 +22,5 @@ fims.log
 ^pkgdown$
 ^assets$
 setup_fims.sh
+^\.positai$
+^\.claude$

.github/workflows/build-metrics.yml

@@ -24,7 +24,7 @@ jobs:
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
       BASELINE_REPOSITORY: ${{ github.event.pull_request.base.repo.full_name || github.repository }}
       CURRENT_COMMIT_SHA: ${{ github.event.pull_request.head.sha || github.sha }}
-      REGRESSION_PCT: 15
+      REGRESSION_PCT: 25
       # Percent threshold for sampled process-tree RSS regressions; sampled peaks are noisier than /usr/bin/time RSS.
       SAMPLED_PROCESS_TREE_RSS_REGRESSION_PCT: 40
     defaults:

.gitignore

@@ -115,3 +115,4 @@ pkgdown/assets/
 
 # CodeQL generated files
 _codeql_detected_source_root
+.positai

CMakeLists.txt

@@ -121,6 +121,25 @@ enable_testing()
 # Creates an interface library fims_test
 add_library(fims_test INTERFACE)
 
+EXECUTE_PROCESS(
+		COMMAND R "--slave" "--no-save" "-e" "cat(system.file('include', package = 'RcppEigen'))"
+		OUTPUT_VARIABLE RCPPEIGEN_INCLUDE_DIR
+		OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+if(NOT RCPPEIGEN_INCLUDE_DIR OR RCPPEIGEN_INCLUDE_DIR STREQUAL "")
+  message(STATUS "RcppEigen not found via R. Attempting to install it...")
+  EXECUTE_PROCESS(
+    COMMAND R "--slave" "--no-save" "-e" "install.packages('RcppEigen', repos='https://cloud.r-project.org'); cat(system.file('include', package = 'RcppEigen'))"
+    OUTPUT_VARIABLE RCPPEIGEN_INCLUDE_DIR
+    OUTPUT_STRIP_TRAILING_WHITESPACE)
+endif()
+
+if(RCPPEIGEN_INCLUDE_DIR AND NOT RCPPEIGEN_INCLUDE_DIR STREQUAL "")
+  message(STATUS "RcppEigen include dir: ${RCPPEIGEN_INCLUDE_DIR}")
+else()
+  message(WARNING "RcppEigen include directory could not be determined. Eigen-dependent tests may fail to compile.")
+endif()
+
 # Add include directories to the fims_test target
 target_include_directories(fims_test
   INTERFACE
@@ -132,8 +151,10 @@ target_include_directories(fims_test
     inst/include/utilities
     ${rcpp_SOURCE_DIR}/inst/include
     ${R_HOME}/include
+    ${RCPPEIGEN_INCLUDE_DIR}
 )
 
+
 # Add compile definition STD_LIB to the fims_test target.
 target_compile_definitions(fims_test
   INTERFACE

NAMESPACE

@@ -33,9 +33,11 @@ export(augment)
 export(clear)
 export(create_default_configurations)
 export(create_default_parameters)
+export(create_edm_embedding)
 export(fit_fims)
 export(get_ages)
 export(get_data)
+export(get_edm_embeddings)
 export(get_end_year)
 export(get_estimates)
 export(get_fit_metrics)
@@ -79,6 +81,7 @@ export(logit)
 export(lognormal)
 export(model_age_comp)
 export(model_age_to_length_conversion)
+export(model_edm_matrix)
 export(model_index)
 export(model_landings)
 export(model_length_comp)

R/Rcpp_exports.R

@@ -8,12 +8,13 @@
 #' [NOAA-FIMS C++ Documentation](https://noaa-fims.github.io/FIMS/doxygen/)
 #'
 #' @name Cpp_classes
-#' @aliases AgeComp BevertonHoltRecruitment CatchAtAge DlnormDistribution DmultinomDistribution DnormDistribution DoubleLogisticSelectivity EWAAGrowth Fleet Index Landings LengthComp LogDevsRecruitmentProcess LogRRecruitmentProcess LogisticMaturity LogisticSelectivity Variable VariableVector Population RealVector SharedInt SharedReal SharedString
-#'
+#' @aliases AgeComp BevertonHoltRecruitment CatchAtAge DelayEmbedding DlnormDistribution DmultinomDistribution DnormDistribution DoubleLogisticSelectivity EWAAGrowth Fleet Index Landings LengthComp LogDevsRecruitmentProcess LogRRecruitmentProcess LogisticMaturity LogisticSelectivity Variable VariableVector Population RealVector SharedInt SharedReal SharedString
+#' 
 #' @details
 #' - [AgeComp](https://noaa-fims.github.io/FIMS/doxygen/classAgeCompDataInterface.html)
 #' - [BevertonHoltRecruitment](https://noaa-fims.github.io/FIMS/doxygen/classBevertonHoltRecruitmentInterface.html)
 #' - [CatchAtAge](https://noaa-fims.github.io/FIMS/doxygen/classCatchAtAgeInterface.html)
+#' - [DelayEmbedding](https://noaa-fims.github.io/FIMS/doxygen/classDelayEmbeddingInterface.html)
 #' - [DlnormDistribution](https://noaa-fims.github.io/FIMS/doxygen/classDlnormDistributionsInterface.html)
 #' - [DmultinomDistribution](https://noaa-fims.github.io/FIMS/doxygen/classDmultinomDistributionsInterface.html)
 #' - [DnormDistribution](https://noaa-fims.github.io/FIMS/doxygen/classDnormDistributionsInterface.html)

R/fimsframe.R

@@ -643,7 +643,7 @@ methods::setMethod(
 #'
 #' @details
 #' The time series is extracted by filtering `get_data(x)` on `type` and
-#' `name` (fleet / survey name). Values encoded as `-999` are treated as the
+#' `fleet` (fleet / survey name). Values encoded as `-999` are treated as the
 #' FIMS missing-data sentinel. When `drop_missing = TRUE`,
 #' `construct_drop_missing()` is called so that rows containing `-999` are
 #' removed from the embedding matrix; when `FALSE`, `construct()` is called
@@ -664,15 +664,15 @@ methods::setMethod(
 #' @param x A [FIMSFrame()] object.
 #' @param series_type A single string giving the `type` value to filter on,
 #'   e.g. `"index"` or `"landings"`.
-#' @param series_name A single string giving the `name` (fleet / survey) value
+#' @param series_name A single string giving the `fleet` (fleet / survey) value
 #'   to filter on, e.g. `"survey1"`.
 #' @param E A positive integer — the embedding dimension (number of lagged
 #'   coordinates per row).
 #' @param tau A positive integer — the time lag between successive coordinates.
 #' @param drop_missing Logical (default `TRUE`). When `TRUE`, rows containing
 #'   the FIMS missing-value sentinel `-999` are dropped from the embedding
 #'   matrix.
-#' @param uncertainty_name A single string giving the `name` column value for
+#' @param uncertainty_name A single string giving the `fleet` column value for
 #'   an uncertainty series in the data slot (e.g. `"survey1_sd"`), or `NULL`
 #'   (default) to skip uncertainty propagation. The uncertainty series must
 #'   have the same length as the value series after filtering on
@@ -720,7 +720,7 @@ create_edm_embedding <- function(
   series_data <- dplyr::filter(
     get_data(x),
     .data[["type"]] == series_type,
-    .data[["name"]] == series_name
+    .data[["fleet"]] == series_name
   ) |>
     dplyr::arrange(.data[["timing"]]) |>
     dplyr::pull(.data[["value"]])
@@ -731,15 +731,15 @@ create_edm_embedding <- function(
     uncertainty_data <- dplyr::filter(
       get_data(x),
       .data[["type"]] == series_type,
-      .data[["name"]] == uncertainty_name
+      .data[["fleet"]] == uncertainty_name
     ) |>
       dplyr::arrange(.data[["timing"]]) |>
       dplyr::pull(.data[["value"]])
 
     if (length(uncertainty_data) == 0) {
       cli::cli_abort(c(
-        "No uncertainty data found for type {.val {series_type}} and name {.val {uncertainty_name}}.",
-        "i" = "Available names: {.val {unique(get_data(x)[[\"name\"]])}}."
+        "No uncertainty data found for type {.val {series_type}} and fleet {.val {uncertainty_name}}.",
+        "i" = "Available fleets: {.val {unique(get_data(x)[[\"fleet\"]])}}."
       ))
     }
     if (length(uncertainty_data) != length(series_data)) {
@@ -813,16 +813,16 @@ create_edm_embedding <- function(
     drop_missing = drop_missing,
     n_rows = as.integer(edm_obj$n_rows),
     n_cols = as.integer(edm_obj$n_cols),
-    embedded_values = edm_obj$embedded_values$toRVector(),
-    target_values = edm_obj$target_values$toRVector()
+    embedded_values = edm_obj$embedded_values$get_values(),
+    target_values = edm_obj$target_values$get_values()
   )
 
   # Append uncertainty fields only when an uncertainty series was provided
   if (!is.null(uncertainty_name)) {
     embedding_result[["embedded_uncertainty"]] <-
-      edm_obj$embedded_uncertainty$toRVector()
+      edm_obj$embedded_uncertainty$get_values()
     embedding_result[["target_uncertainty"]] <-
-      edm_obj$target_uncertainty$toRVector()
+      edm_obj$target_uncertainty$get_values()
   }
 
   # Store in the edm_embeddings slot and return a new FIMSFrame

R/initialize_modules.R

@@ -704,6 +704,8 @@ initialize_fims <- function(parameters, data) {
 
     # TODO(EDM): Detect any standalone EDM likelihood module from the configuration.
 
+    # TODO(EDM): Detect any standalone EDM likelihood module from the configuration.
+
     # Initialize landings module if the data type includes "landings" and
     # if "Landings" exists in the data distribution specification
     if ("landings" %in% fleet_types &&

inst/include/common/fims_math.hpp

@@ -62,6 +62,19 @@ template <class Type>
 inline const Type lgamma(const Type &x) {
   return std::lgamma(x);
 }
+
+/**
+ * @brief The sign function returning -1, 0, or +1.
+ * @param x The value to check the sign of.
+ * @return -1 if x < 0, +1 if x > 0, 0 if x == 0.
+ */
+template <class Type>
+inline const Type sign(const Type &x) {
+  return (x > static_cast<Type>(0))
+             ? static_cast<Type>(1)
+             : ((x < static_cast<Type>(0)) ? static_cast<Type>(-1)
+                                           : static_cast<Type>(0));
+}
 #endif
 
 #ifdef TMB_MODEL
@@ -163,6 +176,17 @@ inline const Type pow(const Type &x, const Type &y) {
   return pow(x, y);
 }
 
+/**
+ * @brief The sign function returning -1, 0, or +1. Uses CppAD::sign under TMB_MODEL.
+ * @param x The value to check the sign of.
+ * @return -1 if x < 0, +1 if x > 0, 0 if x == 0.
+ */
+template <class Type>
+inline const Type sign(const Type &x) {
+  return CppAD::sign(x);
+}
+
+
 /**
  * @brief Computes the natural logarithm of the absolute value of the [gamma
  * function](https://en.wikipedia.org/wiki/Gamma_function) of x for a TMB

inst/include/distributions/distributions.hpp

@@ -11,9 +11,11 @@
 #ifndef FIMS_DISTRIBUTIONS_HPP
 #define FIMS_DISTRIBUTIONS_HPP
 
+#include "functors/beta_lpdf.hpp"
 #include "functors/density_components_base.hpp"
 #include "functors/lognormal_lpdf.hpp"
 #include "functors/multinomial_lpmf.hpp"
+#include "functors/mvnorm_lpdf.hpp"
 #include "functors/normal_lpdf.hpp"
 
 #endif /* FIMS_DISTRIBUTIONS_HPP */

inst/include/distributions/functors/beta_lpdf.hpp

@@ -0,0 +1,108 @@
+/**
+ * @file beta_lpdf.hpp
+ * @brief Implements the BetaLPDF distribution functor used by FIMS to evaluate
+ * observation-level and total log-likelihood contributions under a Beta
+ * distribution model for data, priors, and random effects.
+ * @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 BETA_LPDF
+#define BETA_LPDF
+
+#include <cstddef>
+#include <stdexcept>
+#include <vector>
+
+#include "../../common/def.hpp"
+#include "../../common/fims_math.hpp"
+#include "../../common/fims_vector.hpp"
+#include "density_components_base.hpp"
+
+namespace fims_distributions {
+
+/**
+ * @copybrief beta_lpdf.hpp
+ *
+ * @details This implementation evaluates the Beta log probability density
+ * function for a random variable \f$x \in (0, s)\f$ with shape parameters
+ * \f$\alpha > 0\f$ (`shape1`), \f$\beta > 0\f$ (`shape2`), and upper bound
+ * \f$s > 0\f$ (`scale`, default 1.0):
+ * \f[
+ * \ln f(x \mid \alpha, \beta, s) = (\alpha - 1) \ln(x/s) + (\beta - 1) \ln(1 - x/s)
+ * - \ln \text{B}(\alpha, \beta) - \ln(s)
+ * \f]
+ * where \f$\text{B}(\alpha, \beta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}\f$.
+ *
+ * Under `TMB_MODEL`, observations are passed to `dbeta(x / s, alpha, beta, true) - log(s)`.
+ * Per-observation contributions are stored in `lpdf_vec`; the summed total is returned by
+ * `evaluate()` and stored in `lpdf`.
+ */
+template <typename Type>
+struct BetaLPDF : public DensityComponentBase<Type> {
+  /** @brief Shape parameter 1 (alpha > 0). */
+  fims::Vector<Type> shape1;
+
+  /** @brief Shape parameter 2 (beta > 0). */
+  fims::Vector<Type> shape2;
+
+  /** @brief Scale / upper bound parameter (scale > 0, default 1.0). */
+  fims::Vector<Type> scale;
+
+  /** @brief Constructor. */
+  BetaLPDF() : DensityComponentBase<Type>() {
+    scale.resize(1);
+    scale[0] = static_cast<Type>(1.0);
+  }
+
+  /** @brief Destructor. */
+  virtual ~BetaLPDF() {}
+
+  /**
+   * @brief Evaluates the Beta log probability density function.
+   * @return Summed log probability density value.
+   */
+  virtual const Type evaluate() {
+    size_t n_x = (this->observed_values.size() > 0) ? this->observed_values.size()
+                                                   : this->get_n_x();
+
+    this->lpdf_vec.resize(n_x);
+    std::fill(this->lpdf_vec.begin(), this->lpdf_vec.end(), static_cast<Type>(0));
+    this->lpdf = static_cast<Type>(0);
+
+    for (size_t i = 0; i < n_x; i++) {
+      Type obs = this->get_observed(i);
+      Type a = shape1.get_force_scalar(i);
+      Type b = shape2.get_force_scalar(i);
+      Type s = scale.size() > 0 ? scale.get_force_scalar(i) : static_cast<Type>(1.0);
+
+      Type u = obs / s;
+
+#ifdef TMB_MODEL
+      if (this->input_type == "data") {
+        if (this->get_observed(i) != this->data_observed_values->na_value) {
+          this->lpdf_vec[i] = dbeta(u, a, b, true) - fims_math::log(s);
+        } else {
+          this->lpdf_vec[i] = static_cast<Type>(0);
+        }
+      } else {
+        this->lpdf_vec[i] = dbeta(u, a, b, true) - fims_math::log(s);
+      }
+#else
+      Type log_beta = fims_math::lgamma(a) + fims_math::lgamma(b) -
+                       fims_math::lgamma(a + b);
+      this->lpdf_vec[i] = (a - static_cast<Type>(1.0)) * fims_math::log(u) +
+                          (b - static_cast<Type>(1.0)) * fims_math::log(static_cast<Type>(1.0) - u) -
+                          log_beta - fims_math::log(s);
+#endif
+
+      this->lpdf += this->lpdf_vec[i];
+    }
+
+    return this->lpdf;
+  }
+};
+
+}  // namespace fims_distributions
+
+#endif  // BETA_LPDF