diff --git a/.Rbuildignore b/.Rbuildignore index 2c126f5f9..1e706e80f 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -22,3 +22,5 @@ fims.log ^pkgdown$ ^assets$ setup_fims.sh +^\.positai$ +^\.claude$ diff --git a/.github/workflows/build-metrics.yml b/.github/workflows/build-metrics.yml index 736a45558..ffa778c12 100644 --- a/.github/workflows/build-metrics.yml +++ b/.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: diff --git a/.gitignore b/.gitignore index 0b7eeee0f..480125fe1 100644 --- a/.gitignore +++ b/.gitignore @@ -115,3 +115,4 @@ pkgdown/assets/ # CodeQL generated files _codeql_detected_source_root +.positai diff --git a/CMakeLists.txt b/CMakeLists.txt index 5ac8b2299..bb62ec446 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/NAMESPACE b/NAMESPACE index 806bd38ef..e1c1c2b2d 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/Rcpp_exports.R b/R/Rcpp_exports.R index cf620bd90..133e11ee9 100644 --- a/R/Rcpp_exports.R +++ b/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) diff --git a/R/fimsframe.R b/R/fimsframe.R index ec6f63d6e..82f0f796e 100644 --- a/R/fimsframe.R +++ b/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,7 +664,7 @@ 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). @@ -672,7 +672,7 @@ methods::setMethod( #' @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 diff --git a/R/initialize_modules.R b/R/initialize_modules.R index 08277588b..531dc725d 100644 --- a/R/initialize_modules.R +++ b/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 && diff --git a/inst/include/common/fims_math.hpp b/inst/include/common/fims_math.hpp index 9c38345ac..ab0ab0a00 100644 --- a/inst/include/common/fims_math.hpp +++ b/inst/include/common/fims_math.hpp @@ -62,6 +62,19 @@ template 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 +inline const Type sign(const Type &x) { + return (x > static_cast(0)) + ? static_cast(1) + : ((x < static_cast(0)) ? static_cast(-1) + : static_cast(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 +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 diff --git a/inst/include/distributions/distributions.hpp b/inst/include/distributions/distributions.hpp index dcd12f85d..4a12a6360 100644 --- a/inst/include/distributions/distributions.hpp +++ b/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 */ diff --git a/inst/include/distributions/functors/beta_lpdf.hpp b/inst/include/distributions/functors/beta_lpdf.hpp new file mode 100644 index 000000000..4f9da9011 --- /dev/null +++ b/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 +#include +#include + +#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 +struct BetaLPDF : public DensityComponentBase { + /** @brief Shape parameter 1 (alpha > 0). */ + fims::Vector shape1; + + /** @brief Shape parameter 2 (beta > 0). */ + fims::Vector shape2; + + /** @brief Scale / upper bound parameter (scale > 0, default 1.0). */ + fims::Vector scale; + + /** @brief Constructor. */ + BetaLPDF() : DensityComponentBase() { + scale.resize(1); + scale[0] = static_cast(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(0)); + this->lpdf = static_cast(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(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(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(1.0)) * fims_math::log(u) + + (b - static_cast(1.0)) * fims_math::log(static_cast(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 diff --git a/inst/include/distributions/functors/mvnorm_lpdf.hpp b/inst/include/distributions/functors/mvnorm_lpdf.hpp new file mode 100644 index 000000000..75fdec07f --- /dev/null +++ b/inst/include/distributions/functors/mvnorm_lpdf.hpp @@ -0,0 +1,142 @@ +/** + * @file mvnorm_lpdf.hpp + * @brief Implements the MVNormLPDF distribution functor used by FIMS to + * evaluate observation-level and total log-likelihood contributions under a + * multivariate normal error 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 MVNORM_LPDF +#define MVNORM_LPDF + +#include +#include +#include +#include +#include + +#include "../../common/def.hpp" +#include "../../common/fims_math.hpp" +#include "../../common/fims_vector.hpp" +#include "density_components_base.hpp" + +namespace fims_distributions { + +/** + * @copybrief mvnorm_lpdf.hpp + * + * @details This implementation evaluates the multivariate normal log probability + * density function for a random vector \f$x \in \mathbb{R}^k\f$ with mean vector + * \f$\mu\f$ and covariance matrix \f$\Sigma\f$: + * \f[ + * \ln f(x \mid \mu, \Sigma) = -\frac{k}{2} \ln(2\pi) - \frac{1}{2} \ln|\Sigma| + * -\frac{1}{2} (x - \mu)^\top \Sigma^{-1} (x - \mu) + * \f] + * Under `TMB_MODEL`, it delegates to TMB's `density::MVNORM_t` when available. + * In standalone execution (or general C++ AD), it solves the linear system + * \f$\Sigma \alpha = (x - \mu)\f$ and computes the log determinant using `Eigen::LDLT`. + */ +template +struct MVNormLPDF : public DensityComponentBase { + /** @brief Covariance matrix stored in flat row-major order (k x k). */ + fims::Vector sigma_mat; + + /** @brief Dimension of the multivariate normal vector (k). */ + size_t k_dim = 0; + + /** @brief Constructor. */ + MVNormLPDF() : DensityComponentBase() {} + + /** @brief Destructor. */ + virtual ~MVNormLPDF() {} + + /** + * @brief Evaluates the multivariate normal 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(); + + if (k_dim == 0) { + k_dim = n_x; + } + + if (sigma_mat.size() != k_dim * k_dim) { + throw std::invalid_argument( + "MVNormLPDF: Covariance matrix size (" + + std::to_string(sigma_mat.size()) + ") does not match k_dim^2 (" + + std::to_string(k_dim * k_dim) + ")."); + } + + this->lpdf_vec.resize(1); + this->lpdf_vec[0] = static_cast(0); + this->lpdf = static_cast(0); + +#ifdef TMB_MODEL + // Build Eigen covariance matrix for TMB MVNORM_t + Eigen::Matrix Sigma(k_dim, k_dim); + for (size_t i = 0; i < k_dim; ++i) { + for (size_t j = 0; j < k_dim; ++j) { + Sigma(i, j) = sigma_mat[i * k_dim + j]; + } + } + density::MVNORM_t mvnorm_density(Sigma); + Eigen::Matrix residual(k_dim); + for (size_t i = 0; i < k_dim; ++i) { + Type obs = (this->input_type == "prior" && this->priors.size() == 1) + ? (*(this->priors[0]))[i] + : this->get_observed(i); + Type exp_val = (this->expected_values.size() > 0) + ? this->get_expected(i) + : static_cast(0); + residual(i) = obs - exp_val; + } + // MVNORM_t returns negative log likelihood + this->lpdf_vec[0] = -mvnorm_density(residual); + this->lpdf = this->lpdf_vec[0]; +#else + // Standalone C++ calculation via Eigen::LDLT + Eigen::Matrix Sigma(k_dim, k_dim); + for (size_t i = 0; i < k_dim; ++i) { + for (size_t j = 0; j < k_dim; ++j) { + Sigma(i, j) = sigma_mat[i * k_dim + j]; + } + } + + Eigen::Matrix res(k_dim); + for (size_t i = 0; i < k_dim; ++i) { + Type obs = (this->input_type == "prior" && this->priors.size() == 1) + ? (*(this->priors[0]))[i] + : this->get_observed(i); + Type exp_val = (this->expected_values.size() > 0) + ? this->get_expected(i) + : static_cast(0); + res(i) = obs - exp_val; + } + + Eigen::LDLT> ldlt(Sigma); + Eigen::Matrix alpha = ldlt.solve(res).eval(); + + Type quad_form = static_cast(0); + for (size_t i = 0; i < k_dim; ++i) { + quad_form += res(i) * alpha(i); + } + + // Log-determinant from D vector of LDLT + Type logdet = ldlt.vectorD().array().abs().log().sum(); + + const double log_2pi = std::log(2.0 * M_PI); + this->lpdf_vec[0] = static_cast(-0.5) * + (quad_form + logdet + static_cast(k_dim * log_2pi)); + this->lpdf = this->lpdf_vec[0]; +#endif + + return this->lpdf; + } +}; + +} // namespace fims_distributions + +#endif // MVNORM_LPDF diff --git a/inst/include/edm/edm.hpp b/inst/include/edm/edm.hpp new file mode 100644 index 000000000..5e37d4065 --- /dev/null +++ b/inst/include/edm/edm.hpp @@ -0,0 +1,31 @@ +/** + * @file edm.hpp + * @brief Umbrella header for all EDM (Empirical Dynamic Modelling) components. + * + * @details Including this single header brings in: + * - DelayEmbeddingMatrix, MakeDelayEmbedding, MakeDelayEmbeddingDropMissing + * - SquaredEuclideanDistance, SimplexWeights, SMapWeights, NormalizeWeights + * - Eigen::LDLT (shared linear algebra factorization) + * - ARDKernelElement, BuildCovarianceMatrix, BuildKStarVector (GP kernels) + * - EDMPredictorBase (abstract base for all prediction functors) + * - SimplexProjection (Sugihara & May, 1990) + * - SMapProjection (Sugihara, 1994) + * - GPEdmProjection (Munch et al. 2017; Rogers 2023 GPEDM) + * + * @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_EDM_HPP +#define FIMS_EDM_HPP + +#include "functors/delay_embedding.hpp" +#include "utilities/edm_distance_weights.hpp" +#include "utilities/edm_linear_algebra.hpp" +#include "utilities/edm_kernels.hpp" +#include "functors/edm_predictor_base.hpp" +#include "functors/simplex_projection.hpp" +#include "functors/smap_projection.hpp" +#include "functors/gp_edm_projection.hpp" + +#endif // FIMS_EDM_HPP diff --git a/inst/include/edm/functors/delay_embedding.hpp b/inst/include/edm/functors/delay_embedding.hpp index 887a8897f..912c51f02 100644 --- a/inst/include/edm/functors/delay_embedding.hpp +++ b/inst/include/edm/functors/delay_embedding.hpp @@ -23,7 +23,8 @@ namespace fims_edm { * copying values. This means: * - embedded_values: pointers to the E lagged coordinates per row * [x_t, x_{t-tau}, ..., x_{t-(E-1)tau}] - * - target_values: pointers to the target time-point value x_t for each row + * - target_values: pointers to the forecast target x_{t+h} for each row + * (with default forecast_horizon h=1, this is x_{t+1}) * * Both point into the same underlying input_values vector, so any update to * that vector (e.g., during TMB optimization iterations) is automatically @@ -51,7 +52,15 @@ struct DelayEmbeddingMatrix { std::vector embedded_values; /** * @brief Pointers to the target time-point value for each row. - * target_values[row] == embedded_values[row * n_cols + 0] (i.e., x_t). + * + * With the default forecast_horizon = 1 (one-step-ahead): + * target_values[row] = &series[target_index + 1] (i.e., x_{t+1}) + * + * With forecast_horizon = h: + * target_values[row] = &series[target_index + h] + * + * This matches the standard EDM formulation where the embedding state + * x_t is used to predict x_{t+h} (Sugihara 1994; Esguerra & Munch 2024). */ std::vector target_values; /** @@ -62,8 +71,9 @@ struct DelayEmbeddingMatrix { */ std::vector embedded_uncertainty; /** - * @brief Pointers to the uncertainty at the target time-point for each row. - * target_uncertainty[row] points to sigma_t for that row. + * @brief Pointers to the uncertainty at the forecast target time-point for each row. + * target_uncertainty[row] points to sigma_{target_index + forecast_horizon}, + * i.e., the uncertainty at the predicted time step (mirrors target_values). * Empty when no uncertainty series is provided. */ std::vector target_uncertainty; @@ -88,28 +98,46 @@ struct DelayEmbeddingMatrix { * Rows are ordered by increasing target time. Columns are ordered as * x_t, x_{t - tau}, ..., x_{t - (E - 1) tau}. * + * The target for each row is x_{t + forecast_horizon}. The default + * forecast_horizon = 1 gives the standard one-step-ahead formulation + * (predicting x_{t+1} from the embedding state x_t), matching both + * Sugihara (1994) and Esguerra & Munch (2024, Section 2.1). + * * Both embedded_values and target_values store pointers into @p series so - * that subsequent modifications to @p series are automatically reflected. + * that subsequent modifications to @p series are automatically reflected + * (e.g., during HMS-map EM iterations that update the latent states). * * If @p uncertainty is non-empty (and must be the same length as @p series), * embedded_uncertainty and target_uncertainty are populated with pointers into * @p uncertainty following the same layout. If @p uncertainty is empty, those * fields are left empty. + * + * @throws std::invalid_argument if embedding_dimension or time_lag is 0, + * if the series is too short to accommodate the lag span plus the + * forecast horizon, or if the uncertainty vector has a different length + * from the series. */ template DelayEmbeddingMatrix MakeDelayEmbedding( const fims::Vector& series, size_t embedding_dimension, size_t time_lag, - const fims::Vector& uncertainty = fims::Vector()) { + const fims::Vector& uncertainty = fims::Vector(), + size_t forecast_horizon = 1) { if (embedding_dimension == 0) { throw std::invalid_argument("embedding_dimension must be greater than 0."); } if (time_lag == 0) { throw std::invalid_argument("time_lag must be greater than 0."); } + if (forecast_horizon == 0) { + throw std::invalid_argument("forecast_horizon must be greater than 0."); + } const size_t lag_span = (embedding_dimension - 1) * time_lag; - if (series.size() <= lag_span) { + // Series must be long enough for the lag span plus the forecast horizon. + // Last row has target at index (n - 1 - lag_span + lag_span) + h = n - 1 + h, + // which requires series.size() >= lag_span + forecast_horizon + 1. + if (series.size() < lag_span + forecast_horizon + 1) { throw std::invalid_argument( "series is too short for the requested delay embedding."); } @@ -121,7 +149,10 @@ DelayEmbeddingMatrix MakeDelayEmbedding( } DelayEmbeddingMatrix embedding; - embedding.n_rows = series.size() - lag_span; + // n_rows: the last row has target_index = n_rows - 1 + lag_span and its + // target pointer is at target_index + forecast_horizon. That must be + // < series.size(), so n_rows = series.size() - lag_span - forecast_horizon. + embedding.n_rows = series.size() - lag_span - forecast_horizon; embedding.n_cols = embedding_dimension; embedding.embedded_values.resize(embedding.n_rows * embedding.n_cols); embedding.target_values.resize(embedding.n_rows); @@ -133,15 +164,17 @@ DelayEmbeddingMatrix MakeDelayEmbedding( for (size_t row = 0; row < embedding.n_rows; row++) { const size_t target_index = row + lag_span; - // target_values[row] points to the x_t value for this row - embedding.target_values[row] = &series[target_index]; + // target_values[row] points to x_{target_index + forecast_horizon}, + // i.e., the value the model is trained to predict for this row. + embedding.target_values[row] = &series[target_index + forecast_horizon]; for (size_t col = 0; col < embedding.n_cols; col++) { embedding.embedded_values[row * embedding.n_cols + col] = &series[target_index - col * time_lag]; } if (has_uncertainty) { - embedding.target_uncertainty[row] = &uncertainty[target_index]; + embedding.target_uncertainty[row] = + &uncertainty[target_index + forecast_horizon]; for (size_t col = 0; col < embedding.n_cols; col++) { embedding.embedded_uncertainty[row * embedding.n_cols + col] = &uncertainty[target_index - col * time_lag]; @@ -156,6 +189,9 @@ DelayEmbeddingMatrix MakeDelayEmbedding( * @brief Build a delay embedding matrix and omit rows containing missing * values. * + * Accepts the same @p forecast_horizon parameter as MakeDelayEmbedding() + * and forwards it unchanged. See MakeDelayEmbedding() for full documentation. + * * If @p uncertainty is non-empty (same length as @p series), the uncertainty * fields are propagated for the retained rows using the same pointer-based * approach as MakeDelayEmbedding. @@ -164,9 +200,10 @@ template DelayEmbeddingMatrix MakeDelayEmbeddingDropMissing( const fims::Vector& series, size_t embedding_dimension, size_t time_lag, const Type& missing_value, - const fims::Vector& uncertainty = fims::Vector()) { - DelayEmbeddingMatrix full_embedding = - MakeDelayEmbedding(series, embedding_dimension, time_lag, uncertainty); + const fims::Vector& uncertainty = fims::Vector(), + size_t forecast_horizon = 1) { + DelayEmbeddingMatrix full_embedding = MakeDelayEmbedding( + series, embedding_dimension, time_lag, uncertainty, forecast_horizon); DelayEmbeddingMatrix embedding; embedding.n_cols = full_embedding.n_cols; diff --git a/inst/include/edm/functors/edm_predictor_base.hpp b/inst/include/edm/functors/edm_predictor_base.hpp new file mode 100644 index 000000000..017840ad1 --- /dev/null +++ b/inst/include/edm/functors/edm_predictor_base.hpp @@ -0,0 +1,137 @@ +/** + * @file edm_predictor_base.hpp + * @brief Declares the abstract base class for all EDM prediction functors. + * @details All EDM prediction algorithms (Simplex Projection, S-map, GP-EDM) + * inherit from EDMPredictorBase and override predict_one() and predict(). + * + * Design notes: + * - The library (training manifold) is held as a const pointer into the + * caller-owned DelayEmbeddingMatrix so that TMB optimization updates to the + * underlying series are reflected automatically without rebuilding. + * - EuclideanDistance is a free function in the fims_edm namespace. It + * operates on raw pointers so it can be called inside AD-traced loops + * without relying on STL algorithms that are not available in TMB. + * + * @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_EDM_PREDICTOR_BASE_HPP +#define FIMS_EDM_PREDICTOR_BASE_HPP + +#include +#include +#include + +#include "../../common/fims_vector.hpp" +#include "delay_embedding.hpp" +#include "../utilities/edm_distance_weights.hpp" + +namespace fims_edm { + +// SquaredEuclideanDistance, SimplexWeights, SMapWeights, and NormalizeWeights +// are provided by edm_distance_weights.hpp (included above). + +/** + * @brief Abstract base class for EDM prediction functors. + * + * @details Concrete subclasses must implement: + * - predict_one(): predict a single future value given a query embedding + * coordinate vector. + * - predict(): predict over all rows of a test DelayEmbeddingMatrix and + * populate the `predictions` field. + * + * The caller sets `library` to point at a training DelayEmbeddingMatrix before + * calling either method. The caller is responsible for the lifetime of the + * pointed-to matrix. Passing a null pointer or an incompatible matrix will + * produce a std::runtime_error or std::invalid_argument. + * + * @tparam Type Numeric scalar type (double or TMB AD scalar). + */ +template +struct EDMPredictorBase { + /** + * @brief Pointer to the training manifold (library embedding). + * Not owned by this struct; must outlive all calls to predict_one/predict. + * Must be set by the caller before calling predict_one() or predict(). + */ + const DelayEmbeddingMatrix* library = nullptr; + + /** @brief Embedding dimension E. Must match library->n_cols. */ + size_t embedding_dimension = 0; + + /** @brief Time lag tau (informational; not used by base directly). */ + size_t time_lag = 1; + + /** + * @brief Output predictions. + * Populated by predict(). Element i corresponds to the predicted value for + * test_embedding row i. + */ + std::vector predictions; + + virtual ~EDMPredictorBase() = default; + + /** + * @brief Predict the next value for a single query point. + * + * @param query_point A vector of length embedding_dimension containing the + * lagged coordinates [x_t, x_{t-tau}, ..., x_{t-(E-1)tau}] for the + * query time step. + * @return Predicted value at the next time step. + * @throws std::runtime_error if library is null. + * @throws std::invalid_argument if query_point.size() != embedding_dimension. + */ + virtual Type predict_one(const fims::Vector& query_point) = 0; + + /** + * @brief Predict for all rows of a test embedding. + * + * Iterates over every row of @p test_embedding, extracts the lagged + * coordinate vector, calls predict_one(), and stores results in + * `predictions`. After this call, predictions.size() == test_embedding.n_rows. + * + * @param test_embedding The query manifold. n_cols must equal + * embedding_dimension (and library->n_cols). + * @throws std::runtime_error if library is null. + * @throws std::invalid_argument if dimensions do not match. + */ + virtual void predict(const DelayEmbeddingMatrix& test_embedding) = 0; + + protected: + /** + * @brief Validate that the library pointer is set and dimensions are + * consistent with embedding_dimension. + * @throws std::runtime_error if library is null. + * @throws std::invalid_argument if library->n_cols != embedding_dimension. + */ + void ValidateLibrary() const { + if (library == nullptr) { + throw std::runtime_error( + "EDMPredictorBase: library pointer is null. " + "Set library before calling predict."); + } + if (library->n_cols != embedding_dimension) { + throw std::invalid_argument( + "EDMPredictorBase: library embedding dimension does not match " + "embedding_dimension field."); + } + } + + /** + * @brief Validate that a query point has the correct dimension. + * @param query_point The query coordinate vector. + * @throws std::invalid_argument if size does not match embedding_dimension. + */ + void ValidateQueryPoint(const fims::Vector& query_point) const { + if (query_point.size() != embedding_dimension) { + throw std::invalid_argument( + "EDMPredictorBase: query_point size does not match " + "embedding_dimension."); + } + } +}; + +} // namespace fims_edm + +#endif // FIMS_EDM_PREDICTOR_BASE_HPP diff --git a/inst/include/edm/functors/gp_edm_projection.hpp b/inst/include/edm/functors/gp_edm_projection.hpp new file mode 100644 index 000000000..fa832b778 --- /dev/null +++ b/inst/include/edm/functors/gp_edm_projection.hpp @@ -0,0 +1,555 @@ +/** + * @file gp_edm_projection.hpp + * @brief Implements the GP-EDM prediction algorithm with ARD priors. + * + * @details GP-EDM (Gaussian Process Empirical Dynamic Modelling) predicts the + * next value of a time series using a Gaussian Process regression fitted to + * the delay embedding library. The kernel uses per-dimension inverse + * length-scale parameters (Automatic Relevance Determination, ARD) that are + * estimated by maximising the log posterior via the Rprop gradient optimizer. + * + * ### Mathematical formulation + * + * Given N library rows \f$X \in \mathbb{R}^{N \times E}\f$ and target vector + * \f$y \in \mathbb{R}^N\f$, the covariance matrix is: + * \f[ + * \Sigma_{ij} = \sigma^2 \exp\!\Bigl(-\sum_{d=1}^{E} + * \phi_d (x_{id}-x_{jd})^2\Bigr) + v_e\,\delta_{ij} + * \f] + * + * Hyperparameters optimized by Rprop: + * - \f$\boldsymbol{\phi} = [\phi_1,\ldots,\phi_E]\f$ — ARD inverse + * length-scales; one per embedding dimension (prior: half-Normal). + * - \f$\sigma^2\f$ — signal variance (prior: Beta-shaped on (0, 5)). + * - \f$v_e\f$ — process noise / nugget variance (prior: Beta-shaped on (0, 5)). + * + * The posterior mean prediction for a query point \f$q\f$ is: + * \f[ + * \hat{y}(q) = k^*(q)^\top \Sigma^{-1} y + * \f] + * where \f$k^*(q)_i = \sigma^2 \exp\!\bigl(-\sum_d \phi_d(q_d-x_{id})^2\bigr)\f$. + * + * ### Usage + * ```cpp + * GPEdmProjection gp; + * gp.library = &train_embedding; + * gp.embedding_dimension = E; + * gp.phi = std::vector(E, 0.5); // initial ARD params + * gp.fit(); // optimize phi, sigma2, ve from library + * gp.predict(test_embedding); // fill gp.predictions + * ``` + * + * @references + * Munch, S. B., Poynor, V., and Arriaza, J. L. (2017). Circumventing + * structural uncertainty: a Bayesian perspective on nonlinear forecasting + * for ecology. Ecological Complexity, 32: 134. + * + * Rogers, T. (2023). GPEDM: Gaussian process empirical dynamic modeling. + * https://github.com/tanyarogers/GPEDM (R reference implementation) + * + * @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_EDM_GP_EDM_PROJECTION_HPP +#define FIMS_EDM_GP_EDM_PROJECTION_HPP + +#include +#include +#include +#include +#include +#include + +#include "../../common/fims_math.hpp" +#include "../utilities/edm_kernels.hpp" +#include "../utilities/edm_linear_algebra.hpp" +#include "edm_predictor_base.hpp" + + +namespace fims_edm { + +/** + * @brief GP-EDM prediction functor with per-dimension ARD length-scale priors. + * + * Inherits from EDMPredictorBase. The caller must set `library` and + * `embedding_dimension` before calling fit() or predict_one() / predict(). + * + * Hyperparameters (`phi`, `sigma2`, `ve`) may be set manually or estimated + * from the library by calling fit(). + * + * @tparam Type Numeric scalar type (double or TMB AD scalar). + */ +template +struct GPEdmProjection : public EDMPredictorBase { + // ----------------------------------------------------------------------- + // Hyperparameters (can be set manually or estimated via fit()) + // ----------------------------------------------------------------------- + + /** + * @brief Per-dimension ARD inverse length-scales (length = embedding_dimension). + * + * Each entry \f$\phi_d \geq 0\f$ controls how quickly the kernel decays + * along embedding dimension d: + * - \f$\phi_d \approx 0\f$: dimension d is nearly irrelevant (dropped by ARD). + * - Large \f$\phi_d\f$: strong locality / high nonlinearity along d. + * + * If empty when fit() is called, initialised to 0.1 per dimension. + */ + std::vector phi; + + /** + * @brief Signal variance \f$\sigma^2 > 0\f$. + * Controls the overall amplitude of the covariance function. + */ + double sigma2 = 1.0; + + /** + * @brief Process noise variance \f$v_e > 0\f$ (nugget). + * Added to the diagonal of the covariance matrix for numerical stability + * and to represent observation / process uncertainty. + */ + double ve = 0.1; + + /** + * @brief Maximum Rprop iterations for hyperparameter optimization. + * Matches the default in the GPEDM R package (Rogers 2023). + */ + size_t max_fit_iterations = 200; + + // ----------------------------------------------------------------------- + // Public interface + // ----------------------------------------------------------------------- + + /** + * @brief Estimate hyperparameters from the library via Rprop MAP optimization. + * + * @details Maximises the log posterior: + * \f[ + * \ln p(\boldsymbol{\phi}, \sigma^2, v_e \mid y) = + * -\tfrac{1}{2} y^\top \Sigma^{-1} y - \tfrac{1}{2} \ln|\Sigma| + * + \ln \pi(\boldsymbol{\phi}) + \ln \pi(\sigma^2) + \ln \pi(v_e) + * \f] + * using the Rprop (Resilient Backpropagation) optimizer with adaptive + * per-parameter step sizes, matching the GPEDM R package (Rogers 2023). + * + * Parameters are optimised in transformed (unconstrained) space: + * - \f$\phi_d\f$: \f$\ln \phi_d\f$ (log transform) + * - \f$\sigma^2\f$: logit on (0.0001, 5) + * - \f$v_e\f$: logit on (0.0001, 5) + * + * After calling fit(), `phi`, `sigma2`, and `ve` contain the MAP estimates. + * + * @throws std::runtime_error if library is null or N < 2. + * @throws std::invalid_argument if embedding_dimension mismatches library. + */ + void fit() { + this->ValidateLibrary(); + const size_t N = this->library->n_rows; + const size_t E = this->embedding_dimension; + + if (N < 2) { + throw std::runtime_error( + "GPEdmProjection::fit: library must have at least 2 rows."); + } + + // --- Initialize phi if empty --- + if (phi.empty()) { + phi.assign(E, 0.1); + } + if (phi.size() != E) { + throw std::invalid_argument( + "GPEdmProjection::fit: phi.size() must equal embedding_dimension."); + } + + // --- Bounds for sigma2 and ve (same as GPEDM R package) --- + const double kVeMin = 0.0001, kVeMax = 4.9999; + const double kS2Min = 0.0001, kS2Max = 4.9999; + + // --- Build initial parameter vector in transformed space --- + // parst = [log(phi[0]), ..., log(phi[E-1]), logit(ve), logit(sigma2)] + const size_t np = E + 2; + std::vector parst(np); + for (size_t d = 0; d < E; ++d) { + parst[d] = std::log(phi[d] > 1e-12 ? phi[d] : 1e-12); + } + parst[E] = logit_transform(ve, kVeMin, kVeMax); + parst[E + 1] = logit_transform(sigma2, kS2Min, kS2Max); + + // --- Extract library targets into a plain double vector --- + std::vector y(N); + for (size_t i = 0; i < N; ++i) { + y[i] = static_cast(*(this->library->target_values[i])); + } + + // --- Pre-compute per-dimension squared-distance matrices (N x N) --- + // D[d][i*N + j] = (x_id - x_jd)^2 + std::vector> D(E, std::vector(N * N, 0.0)); + for (size_t d = 0; d < E; ++d) { + for (size_t i = 0; i < N; ++i) { + double xi = + static_cast(*(this->library->embedded_values[i * E + d])); + for (size_t j = i; j < N; ++j) { + double xj = + static_cast(*(this->library->embedded_values[j * E + d])); + double diff = xi - xj; + double sq = diff * diff; + D[d][i * N + j] = sq; + D[d][j * N + i] = sq; + } + } + } + + // --- Rprop optimizer --- + // Adaptive per-parameter step sizes, sign-based descent. + const double kEtaPlus = 0.2; // step growth (factor 1+eta) + const double kEtaMinus = -0.5; // step shrink (factor 1+eta) + const double kDeltaMin = 1e-6; + const double kDeltaMax = 50.0; + + std::vector delta(np, 0.1); + std::vector grad_prev(np, 0.0); + double nllpost_prev = std::numeric_limits::max(); + + // Compute initial gradient + LogPosteriorResult res = compute_log_posterior_grad( + parst, y, D, N, E, kVeMin, kVeMax, kS2Min, kS2Max); + + for (size_t iter = 0; iter < max_fit_iterations; ++iter) { + // Convergence check (gradient norm and relative change in objective) + double gnorm = 0.0; + for (size_t k = 0; k < np; ++k) gnorm += res.grad[k] * res.grad[k]; + gnorm = std::sqrt(gnorm); + double df = std::abs(res.nllpost / (nllpost_prev + 1e-15) - 1.0); + if (gnorm < 1e-4 && df < 1e-7) break; + nllpost_prev = res.nllpost; + + // --- Rprop step --- + std::vector parst_new(np); + for (size_t k = 0; k < np; ++k) { + parst_new[k] = parst[k] - fims_math::sign(res.grad[k]) * delta[k]; + } + + + LogPosteriorResult res_new = compute_log_posterior_grad( + parst_new, y, D, N, E, kVeMin, kVeMax, kS2Min, kS2Max); + + // Update adaptive step sizes + for (size_t k = 0; k < np; ++k) { + double gc = res.grad[k] * res_new.grad[k]; + double tdelta = + delta[k] * (1.0 + (gc > 0.0 ? kEtaPlus : (gc < 0.0 ? kEtaMinus : 0.0))); + delta[k] = tdelta < kDeltaMin ? kDeltaMin + : (tdelta > kDeltaMax ? kDeltaMax : tdelta); + } + + parst = parst_new; + res = res_new; + } + + // --- Untransform optimized parameters --- + for (size_t d = 0; d < E; ++d) { + phi[d] = std::exp(parst[d]); + } + ve = logit_inverse(parst[E], kVeMin, kVeMax); + sigma2 = logit_inverse(parst[E + 1], kS2Min, kS2Max); + } + + /** + * @brief Predict the next value for a single query embedding point. + * + * @details Computes the GP posterior mean: + * \f[ + * \hat{y}(q) = k^*(q)^\top \Sigma^{-1} y + * \f] + * using the current values of `phi`, `sigma2`, and `ve`. If fit() has + * not been called, the default hyperparameters are used directly. + * + * @param query_point Length-E vector of lagged coordinates for the query. + * @return Predicted value at the next time step. + * @throws std::runtime_error if library is null or has fewer than 2 rows. + * @throws std::invalid_argument if query_point.size() != embedding_dimension. + */ + Type predict_one(const fims::Vector& query_point) override { + this->ValidateLibrary(); + this->ValidateQueryPoint(query_point); + + const size_t N = this->library->n_rows; + const size_t E = this->embedding_dimension; + + if (N < 2) { + throw std::runtime_error( + "GPEdmProjection::predict_one: library must have at least 2 rows."); + } + if (phi.size() != E) { + throw std::invalid_argument( + "GPEdmProjection::predict_one: phi.size() must equal " + "embedding_dimension. Call fit() first or set phi manually."); + } + + // --- Build N×N covariance matrix Sigma = K + ve*I --- + std::vector Sigma; + BuildCovarianceMatrix(this->library->embedded_values, N, E, + phi, sigma2, ve, Sigma); + + // --- Build target vector y --- + std::vector y(N); + for (size_t i = 0; i < N; ++i) { + y[i] = *(this->library->target_values[i]); + } + + // --- Solve Sigma * alpha = y (alpha overwrites y) via Cholesky (Eigen::LDLT) --- + Eigen::Matrix Sigma_mat(N, N); + for (size_t i = 0; i < N; ++i) { + for (size_t j = 0; j < N; ++j) { + Sigma_mat(i, j) = Sigma[i * N + j]; + } + } + Eigen::Matrix y_mat(N); + for (size_t i = 0; i < N; ++i) { + y_mat(i) = y[i]; + } + Eigen::LDLT> ldlt(Sigma_mat); + Eigen::Matrix alpha = ldlt.solve(y_mat).eval(); + for (size_t i = 0; i < N; ++i) { + y[i] = alpha(i); + } + // y now holds alpha = Sigma^{-1} y_targets + + + // --- Build k_star: covariance between query and every library row --- + std::vector k_star; + BuildKStarVector(query_point.data(), this->library->embedded_values, + N, E, phi, sigma2, k_star); + + // --- Posterior mean: y_hat = k_star^T alpha --- + Type y_hat = Type(0); + for (size_t i = 0; i < N; ++i) { + y_hat += k_star[i] * y[i]; + } + return y_hat; + } + + /** + * @brief Predict for every row of a test embedding matrix. + * + * @details Iterates over @p test_embedding rows, extracts each query vector, + * calls predict_one(), and stores results in `predictions`. + * + * @param test_embedding Query manifold (n_cols must equal embedding_dimension). + * @throws std::runtime_error if library is null or has fewer than 2 rows. + * @throws std::invalid_argument if test_embedding.n_cols != embedding_dimension. + */ + void predict(const DelayEmbeddingMatrix& test_embedding) override { + this->ValidateLibrary(); + if (test_embedding.n_cols != this->embedding_dimension) { + throw std::invalid_argument( + "GPEdmProjection::predict: test_embedding.n_cols does not match " + "embedding_dimension."); + } + + const size_t n_test = test_embedding.n_rows; + this->predictions.resize(n_test); + + fims::Vector query(this->embedding_dimension); + for (size_t row = 0; row < n_test; ++row) { + for (size_t col = 0; col < this->embedding_dimension; ++col) { + query[col] = test_embedding.at(row, col); + } + this->predictions[row] = predict_one(query); + } + } + + private: + // ----------------------------------------------------------------------- + // Rprop helpers (all in double arithmetic, not in the AD trace) + // ----------------------------------------------------------------------- + + /** @brief Result of one log-posterior evaluation. */ + struct LogPosteriorResult { + double nllpost; ///< Negative log posterior value. + std::vector grad; ///< Gradient of nllpost w.r.t. parst. + }; + + /** @brief Logit transform: maps x in (lo, hi) to the real line. */ + static double logit_transform(double x, double lo, double hi) { + double p = (x - lo) / (hi - lo); + p = p < 1e-9 ? 1e-9 : (p > 1.0 - 1e-9 ? 1.0 - 1e-9 : p); + return std::log(p / (1.0 - p)); + } + + /** @brief Inverse logit: maps real line back to (lo, hi). */ + static double logit_inverse(double t, double lo, double hi) { + return (hi - lo) / (1.0 + fims_math::exp(-t)) + lo; + } + + + /** + * @brief Evaluate the negative log posterior and its gradient. + * + * @details This runs entirely in double arithmetic (not on the AD tape). + * Uses the analytic gradient derived from the GP marginal likelihood: + * \f{align*}{ + * \ln p &= -\tfrac{1}{2} y^\top \Sigma^{-1} y + * -\tfrac{1}{2} \ln|\Sigma| + \ln\pi(\boldsymbol{\phi}, + * \sigma^2, v_e) \\ + * \nabla_{\phi_d} \ln p &= \tfrac{1}{2} \operatorname{tr} + * \bigl((\alpha\alpha^\top - \Sigma^{-1}) \partial_{\phi_d}\Sigma\bigr)\\ + * \partial_{\phi_d}\Sigma_{ij} &= -D_{d,ij}\, K_{ij} + * \f} + * where \f$\alpha = \Sigma^{-1} y\f$ and \f$D_{d,ij} = (x_{id}-x_{jd})^2\f$. + * + * @param parst Transformed parameters (log phi; logit ve; logit sigma2). + * @param y Target vector (length N, plain double). + * @param D Per-dimension squared-distance matrices (E × N*N). + * @param N, E Library dimensions. + * @param ve_min, ve_max, s2_min, s2_max Parameter bounds. + */ + LogPosteriorResult compute_log_posterior_grad( + const std::vector& parst, const std::vector& y, + const std::vector>& D, size_t N, size_t E, + double ve_min, double ve_max, double s2_min, double s2_max) const { + const size_t np = E + 2; + + // --- Untransform parameters --- + std::vector phi_cur(E); + for (size_t d = 0; d < E; ++d) phi_cur[d] = std::exp(parst[d]); + double ve_cur = logit_inverse(parst[E], ve_min, ve_max); + double sigma2_cur = logit_inverse(parst[E + 1], s2_min, s2_max); + + // --- Jacobian factors (chain rule through transforms) --- + // d(nllpost)/d(parst[d]) = d(nllpost)/d(phi[d]) * d(phi[d])/d(parst[d]) + // d(phi[d])/d(log phi[d]) = phi[d] + // d(ve)/d(logit ve) = (ve - ve_min)(1 - (ve - ve_min)/(ve_max - ve_min)) + std::vector dpars(np, 0.0); + for (size_t d = 0; d < E; ++d) dpars[d] = phi_cur[d]; + dpars[E] = (ve_cur - ve_min) * (1.0 - (ve_cur - ve_min) / (ve_max - ve_min)); + dpars[E + 1] = (sigma2_cur - s2_min) * (1.0 - (sigma2_cur - s2_min) / (s2_max - s2_min)); + + // --- Build covariance matrix Sigma = K + ve*I (double version) --- + std::vector Sigma_d(N * N, 0.0); + std::vector K0(N * N, 0.0); // K without nugget, for gradient + for (size_t i = 0; i < N; ++i) { + for (size_t j = i; j < N; ++j) { + double sq = 0.0; + for (size_t d = 0; d < E; ++d) sq += phi_cur[d] * D[d][i * N + j]; + double k_ij = sigma2_cur * std::exp(-sq); + K0[i * N + j] = k_ij; + K0[j * N + i] = k_ij; + Sigma_d[i * N + j] = k_ij; + Sigma_d[j * N + i] = k_ij; + } + Sigma_d[i * N + i] += ve_cur; + } + + // --- Solve Sigma * alpha = y and compute logdet & iKVs via Eigen::LDLT --- + Eigen::MatrixXd Sigma_mat(N, N); + for (size_t i = 0; i < N; ++i) { + for (size_t j = 0; j < N; ++j) { + Sigma_mat(i, j) = Sigma_d[i * N + j]; + } + } + Eigen::VectorXd y_eig(N); + for (size_t i = 0; i < N; ++i) { + y_eig(i) = y[i]; + } + + Eigen::LDLT ldlt(Sigma_mat); + Eigen::VectorXd alpha_eig = ldlt.solve(y_eig).eval(); + std::vector alpha(N); + for (size_t i = 0; i < N; ++i) { + alpha[i] = alpha_eig(i); + } + + // Inverse matrix iKVs = Sigma^{-1} + Eigen::MatrixXd I_mat = Eigen::MatrixXd::Identity(N, N); + Eigen::MatrixXd iKVs_mat = ldlt.solve(I_mat).eval(); + std::vector iKVs(N * N, 0.0); + for (size_t i = 0; i < N; ++i) { + for (size_t j = 0; j < N; ++j) { + iKVs[i * N + j] = iKVs_mat(i, j); + } + } + + // Log determinant from LDLT vectorD + double logdet = ldlt.vectorD().array().abs().log().sum(); + + double ytAlpha = 0.0; + for (size_t i = 0; i < N; ++i) ytAlpha += y[i] * alpha[i]; + double like = -0.5 * ytAlpha - 0.5 * logdet; + + + // --- Priors (matching GPEDM R package: Rogers 2023, getpriors()) --- + // phi: half-Normal prior + // lam_phi = (2^(modeprior-1))^2 * pi/2 + // Default modeprior=0: lam = (2^(-1))^2 * pi/2 = 0.25 * pi/2 = pi/8 + const double kLamPhi = 0.25 * M_PI / 2.0; // = pi/8 + double lp_phi = 0.0; + std::vector dlp_phi(E, 0.0); + for (size_t d = 0; d < E; ++d) { + lp_phi += -0.5 * phi_cur[d] * phi_cur[d] / kLamPhi; + dlp_phi[d] = -phi_cur[d] / kLamPhi; + } + // sigma2 and ve: Beta-shaped (a=2, b=2) on (0, max) + // lp = (a-1)*log(x/xmax) + (b-1)*log(1 - x/xmax) with a=b=2 → (a-1)=1 + double lp_ve = std::log(ve_cur / ve_max) + std::log(1.0 - ve_cur / ve_max); + double lp_s2 = std::log(sigma2_cur / s2_max) + std::log(1.0 - sigma2_cur / s2_max); + // dlp = (a-1)/x - (b-1)/(xmax-x) with a=b=2 → (a-1)=1 + double dlp_ve = 1.0 / ve_cur - 1.0 / (ve_max - ve_cur); + double dlp_s2 = 1.0 / sigma2_cur - 1.0 / (s2_max - sigma2_cur); + + double lp = lp_phi + lp_ve + lp_s2; + double nllpost = -(like + lp); + + // --- Gradient of log posterior --- + // vQ[i,j] = 0.5 * (alpha[i]*alpha[j] - iKVs[i,j]) + // dl[d] = 0.5 * sum_{i,j} vQ[i,j] * dSigma[i,j]/dphi[d] + // = 0.5 * sum_{i,j} vQ[i,j] * (-D[d][i,j] * K0[i,j]) + std::vector grad_unconstrained(np, 0.0); + + for (size_t d = 0; d < E; ++d) { + double dl_d = 0.0; + for (size_t i = 0; i < N; ++i) { + for (size_t j = 0; j < N; ++j) { + double vQ_ij = 0.5 * (alpha[i] * alpha[j] - iKVs[i * N + j]); + dl_d += vQ_ij * (-D[d][i * N + j] * K0[i * N + j]); + } + } + grad_unconstrained[d] = dl_d + dlp_phi[d]; + } + + // ve gradient: dSigma/dve = I + { + double dl_ve = 0.0; + for (size_t i = 0; i < N; ++i) { + dl_ve += 0.5 * (alpha[i] * alpha[i] - iKVs[i * N + i]); + } + grad_unconstrained[E] = dl_ve + dlp_ve; + } + + // sigma2 gradient: dSigma/dsigma2 = K0 / sigma2 + { + double dl_s2 = 0.0; + for (size_t i = 0; i < N; ++i) { + for (size_t j = 0; j < N; ++j) { + double vQ_ij = 0.5 * (alpha[i] * alpha[j] - iKVs[i * N + j]); + dl_s2 += vQ_ij * (K0[i * N + j] / (sigma2_cur + 1e-300)); + } + } + grad_unconstrained[E + 1] = dl_s2 + dlp_s2; + } + + // Chain-rule through parameter transforms: grad_parst = -J * dpars + LogPosteriorResult result; + result.nllpost = nllpost; + result.grad.resize(np); + for (size_t k = 0; k < np; ++k) { + result.grad[k] = -grad_unconstrained[k] * dpars[k]; + } + return result; + } +}; + +} // namespace fims_edm + +#endif // FIMS_EDM_GP_EDM_PROJECTION_HPP diff --git a/inst/include/edm/functors/simplex_projection.hpp b/inst/include/edm/functors/simplex_projection.hpp new file mode 100644 index 000000000..49ec78e83 --- /dev/null +++ b/inst/include/edm/functors/simplex_projection.hpp @@ -0,0 +1,194 @@ +/** + * @file simplex_projection.hpp + * @brief Implements the Simplex Projection EDM prediction algorithm. + * + * @details Simplex Projection (Sugihara & May, 1990) predicts the next value + * of a time series by finding the E+1 nearest neighbors of the query point in + * the library manifold and returning an exponentially-weighted average of + * their one-step-ahead targets. + * + * ### Algorithm + * Given a query embedding coordinate `q` of dimension E: + * + * 1. Compute the squared Euclidean distance d_i from q to every library row i. + * 2. Select the k = E+1 nearest neighbors (hand-rolled O(n*k) partial + * selection to remain compatible with TMB's CppAD tracer, which cannot + * trace std::partial_sort or std::nth_element). + * 3. Compute exponential weights: + * \f[ + * w_i = \exp\!\Bigl(-\frac{d_i}{d_{\min} + \varepsilon}\Bigr) + * \f] + * where d_min is the smallest neighbor distance and eps guards against + * division by zero when the query is an exact match. + * 4. Return the normalized weighted average of the library target values + * for those k neighbors: + * \f[ + * \hat{x}_{t+1} = \frac{\sum_i w_i \cdot x^*_i}{\sum_i w_i} + * \f] + * + * @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_EDM_SIMPLEX_PROJECTION_HPP +#define FIMS_EDM_SIMPLEX_PROJECTION_HPP + +#include +#include +#include +#include + +#include "../../common/fims_math.hpp" +#include "../utilities/edm_distance_weights.hpp" +#include "edm_predictor_base.hpp" + +namespace fims_edm { + +/** + * @brief Simplex Projection EDM predictor. + * + * @details Inherits from EDMPredictorBase. Set `library` and + * `embedding_dimension` before calling predict_one() or predict(). + * + * `n_neighbors` defaults to embedding_dimension + 1 when left at 0 (the + * theoretically motivated choice from Sugihara & May 1990). The caller may + * override it for smoother (larger k) or sharper (smaller k) predictions. + * + * @tparam Type Numeric scalar type (double or a TMB AD scalar). + */ +template +struct SimplexProjection : public EDMPredictorBase { + /** + * @brief Number of nearest neighbors. + * When 0 (default), predict_one() uses embedding_dimension + 1. + * Override to a positive value to fix k regardless of E. + */ + size_t n_neighbors = 0; + + /** + * @brief Small constant added to d_min when computing exponential weights. + * Prevents NaN when d_min == 0 (exact query match in the library). + */ + double weight_epsilon = 1e-12; + + SimplexProjection() = default; + virtual ~SimplexProjection() = default; + + /** + * @brief Predict the next value for a single query embedding point. + * + * @param query_point Lagged coordinate vector of length embedding_dimension + * representing the state at time t: + * [x_t, x_{t-tau}, ..., x_{t-(E-1)*tau}]. + * @return Predicted value \f$\hat{x}_{t+1}\f$. + * @throws std::runtime_error if library is null. + * @throws std::invalid_argument if query_point.size() != embedding_dimension. + * @throws std::runtime_error if the library has fewer rows than n_neighbors. + */ + virtual Type predict_one(const fims::Vector& query_point) override { + this->ValidateLibrary(); + this->ValidateQueryPoint(query_point); + + const size_t k = + (n_neighbors == 0) ? (this->embedding_dimension + 1) : n_neighbors; + const size_t n_lib = this->library->n_rows; + const size_t E = this->embedding_dimension; + + if (n_lib < k) { + throw std::runtime_error( + "SimplexProjection: library has fewer rows than n_neighbors (k). " + "Increase the library size or decrease n_neighbors."); + } + + // ----------------------------------------------------------------------- + // Step 1 & 2: Compute squared distances and track the k nearest neighbors + // with a hand-rolled O(n*k) partial selection. Using raw pointer indexing + // into embedded_values so the loop is CppAD-traceable. + // + // nn_indices[m] = library row index of the m-th nearest-neighbor slot. + // nn_sq_dist[m] = its squared distance (initialized to +inf). + // ----------------------------------------------------------------------- + const double kInf = std::numeric_limits::max(); + std::vector nn_indices(k, 0); + std::vector nn_sq_dist(k, Type(kInf)); + + for (size_t i = 0; i < n_lib; ++i) { + // Squared distance from query_point to library row i. + Type d = Type(0); + for (size_t j = 0; j < E; ++j) { + Type diff = query_point[j] - *(this->library->embedded_values[i * E + j]); + d += diff * diff; + } + + // Find the slot with the current worst (largest) distance. + size_t worst = 0; + for (size_t m = 1; m < k; ++m) { + if (nn_sq_dist[m] > nn_sq_dist[worst]) { + worst = m; + } + } + + // Replace the worst slot if this neighbor is closer. + if (d < nn_sq_dist[worst]) { + nn_sq_dist[worst] = d; + nn_indices[worst] = i; + } + } + + // ----------------------------------------------------------------------- + // Step 3: Compute Simplex exponential weights using the shared utility. + // w_i = exp(-d_i / (d_min + eps)), then normalize so they sum to 1. + // ----------------------------------------------------------------------- + std::vector weights; + SimplexWeights(nn_sq_dist, weights, weight_epsilon); + NormalizeWeights(weights); + + // ----------------------------------------------------------------------- + // Step 4: Normalized weighted average of target values. + // ----------------------------------------------------------------------- + Type prediction = Type(0); + for (size_t m = 0; m < k; ++m) { + prediction += weights[m] * (*this->library->target_values[nn_indices[m]]); + } + return prediction; + } + + /** + * @brief Predict for all rows of a test embedding matrix. + * + * Iterates over every row of @p test_embedding, builds the query coordinate + * vector, calls predict_one(), and stores the results in `predictions`. + * After this call: predictions.size() == test_embedding.n_rows. + * + * @param test_embedding The query manifold. n_cols must equal + * embedding_dimension. + * @throws std::runtime_error / std::invalid_argument forwarded from + * predict_one() and ValidateLibrary(). + */ + virtual void predict( + const DelayEmbeddingMatrix& test_embedding) override { + this->ValidateLibrary(); + + const size_t E = this->embedding_dimension; + if (test_embedding.n_cols != E) { + throw std::invalid_argument( + "SimplexProjection::predict: test embedding n_cols does not match " + "embedding_dimension."); + } + + const size_t n_test = test_embedding.n_rows; + this->predictions.resize(n_test); + + fims::Vector query_point(E); + for (size_t row = 0; row < n_test; ++row) { + for (size_t col = 0; col < E; ++col) { + query_point[col] = test_embedding.at(row, col); + } + this->predictions[row] = predict_one(query_point); + } + } +}; + +} // namespace fims_edm + +#endif // FIMS_EDM_SIMPLEX_PROJECTION_HPP diff --git a/inst/include/edm/functors/smap_projection.hpp b/inst/include/edm/functors/smap_projection.hpp new file mode 100644 index 000000000..796040efb --- /dev/null +++ b/inst/include/edm/functors/smap_projection.hpp @@ -0,0 +1,266 @@ +/** + * @file smap_projection.hpp + * @brief Implements the S-Map (Sequential Locally Weighted Global Linear Maps) + * EDM prediction algorithm. + * + * @details S-Map (Sugihara 1994) predicts the next value of a time series by + * fitting a **locally weighted linear regression** using ALL library rows. + * Unlike Simplex Projection (which uses only E+1 nearest neighbors), + * S-Map uses every library point but down-weights distant ones via a + * kernel function controlled by the nonlinearity parameter theta (θ). + * + * Two kernel functions are supported (see SMapKernel): + * - SMapKernel::kExponential (default): w = exp(-θ · d / d̄) [rEDM style] + * - SMapKernel::kGaussian: w = exp(-θ² · (d/D)²) [Esguerra & Munch 2024] + * Per Munch (pers. comm.) the performance difference is rarely large; + * use kGaussian to avoid a sqrt() inside the AD tape. + * + * ### Algorithm + * Given a query embedding coordinate q of dimension E: + * + * 1. Compute squared Euclidean distance d_i from q to every library row i. + * 2. Compute S-Map weights using SMapWeights() from edm_distance_weights.hpp: + * \f[ + * w_i = \exp\!\Bigl(-\theta \cdot \frac{d_i}{\bar{d}}\Bigr) + * \f] + * where \f$\bar{d}\f$ is the mean squared distance across all library rows + * and \f$\theta \geq 0\f$ controls locality: + * - \f$\theta = 0\f$: all weights equal → global linear model. + * - Large \f$\theta\f$: only nearby rows contribute → highly local. + * 3. Assemble the weighted normal equations for coefficient vector + * \f$\boldsymbol{\beta}\f$ of size E+1 (intercept + E slopes): + * \f[ + * \mathbf{A} = \mathbf{X}^{\top} \mathbf{W} \mathbf{X},\quad + * \mathbf{b} = \mathbf{X}^{\top} \mathbf{W} \mathbf{y} + * \f] + * where row i of \f$\mathbf{X}\f$ is \f$[1,\, x_{i1},\ldots, x_{iE}]\f$, + * \f$\mathbf{W}\f$ is diagonal with entries \f$w_i\f$, and + * \f$y_i\f$ is the library target value. + * 4. Solve \f$\mathbf{A}\boldsymbol{\beta} = \mathbf{b}\f$ via Gaussian + * elimination with partial pivoting (in-place, CppAD-compatible). + * 5. Return \f$\hat{x}_{t+1} = [1, q_1, \ldots, q_E]\cdot\boldsymbol{\beta}\f$. + * + * ### Relationship to Simplex Projection + * | Feature | Simplex | S-Map | + * |---|---|---| + * | Library rows used | E+1 nearest | ALL | + * | Prediction model | weighted average | locally weighted linear regression | + * | Key parameter | n_neighbors k | nonlinearity θ | + * | Weighting (default) | exp(-d/d_min) | exp(-θ·d/d̄) | + * | Weighting (Gaussian) | — | exp(-θ²·(d/D)²) | + * + * @see edm_distance_weights.hpp for SMapWeights() and NormalizeWeights(). + * + * @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_EDM_SMAP_PROJECTION_HPP +#define FIMS_EDM_SMAP_PROJECTION_HPP + +#include +#include +#include +#include + +#include "../utilities/edm_distance_weights.hpp" +#include "../utilities/edm_linear_algebra.hpp" +#include "edm_predictor_base.hpp" + +namespace fims_edm { + +// --------------------------------------------------------------------------- +// SMapProjection +// --------------------------------------------------------------------------- + +/** + * @brief S-Map EDM predictor (Sugihara 1994). + * + * @details Inherits from EDMPredictorBase. Set `library` and + * `embedding_dimension` before calling predict_one() or predict(). + * + * @tparam Type Numeric scalar type (double or a TMB AD scalar). + */ +template +struct SMapProjection : public EDMPredictorBase { + /** + * @brief Nonlinearity parameter θ ≥ 0. + * + * Controls how strongly the model localizes around the query point. + * - θ = 0 → global linear regression (all weights equal 1). + * - θ > 0 → down-weight distant library rows via the chosen kernel. + * Typical values: 0, 0.5, 1, 2, 4, 8 (scan via cross-validation). + */ + double theta = 1.0; + + /** + * @brief Weighting kernel used to compute S-Map locality weights. + * + * Choices: + * - SMapKernel::kExponential (default): w = exp(-θ·d/d̄) [classic rEDM] + * - SMapKernel::kGaussian: w = exp(-θ²·(d/D)²) [Esguerra & Munch 2024] + * + * Both kernels reduce to uniform weights when θ = 0 (global OLS). + * Per Munch (pers. comm.) the performance difference is rarely large; + * kGaussian avoids a sqrt() call inside the CppAD tape. + */ + SMapKernel kernel = SMapKernel::kExponential; + + /** + * @brief Small constant added to the mean distance to prevent division + * by zero when all library rows are identical to the query. + */ + double epsilon = 1e-12; + + SMapProjection() = default; + virtual ~SMapProjection() = default; + + /** + * @brief Predict the next value for a single query embedding point. + * + * @details Fits a locally weighted linear regression using all library rows, + * with weights w_i = exp(-θ * d_i / d̄). The regression includes an + * intercept term, so the coefficient vector β has length E+1. + * + * @param query_point Lagged coordinate vector of length embedding_dimension + * representing the state at time t: + * [x_t, x_{t-tau}, ..., x_{t-(E-1)*tau}]. + * @return Predicted value \f$\hat{x}_{t+1}\f$. + * @throws std::runtime_error if library is null. + * @throws std::invalid_argument if query_point.size() != embedding_dimension. + * @throws std::runtime_error if the weighted design matrix is singular. + */ + virtual Type predict_one(const fims::Vector& query_point) override { + this->ValidateLibrary(); + this->ValidateQueryPoint(query_point); + + const size_t n_lib = this->library->n_rows; + const size_t E = this->embedding_dimension; + + if (n_lib < E + 1) { + throw std::runtime_error( + "SMapProjection: library must have at least E+1 rows to fit a " + "linear model with E predictors and an intercept."); + } + + // ----------------------------------------------------------------------- + // Step 1: Compute squared distance from query to ALL library rows. + // Raw pointer indexing keeps the loop CppAD-traceable. + // ----------------------------------------------------------------------- + std::vector sq_dists(n_lib, Type(0)); + for (size_t i = 0; i < n_lib; ++i) { + for (size_t j = 0; j < E; ++j) { + Type diff = query_point[j] - + *(this->library->embedded_values[i * E + j]); + sq_dists[i] += diff * diff; + } + } + + // ----------------------------------------------------------------------- + // Step 2: Compute S-Map weights using the user-selected kernel. + // kExponential: w_i = exp(-theta * d_i / d_mean) + // kGaussian: w_i = exp(-theta^2 * (d_i / D)^2) + // ----------------------------------------------------------------------- + std::vector weights; + SMapWeights(sq_dists, weights, theta, kernel, epsilon); + + // ----------------------------------------------------------------------- + // Step 3: Assemble normal equations A β = b. + // + // Design matrix row i: x̃_i = [1, lib[i][0], ..., lib[i][E-1]] (length E+1) + // A[k][j] += w_i * x̃_i[k] * x̃_i[j] + // b[k] += w_i * x̃_i[k] * y_i + // + // A is stored flat, row-major, size (E+1)*(E+1). + // ----------------------------------------------------------------------- + const size_t p = E + 1; // number of parameters (intercept + E slopes) + std::vector A(p * p, Type(0)); + std::vector b_vec(p, Type(0)); + + for (size_t i = 0; i < n_lib; ++i) { + const Type w = weights[i]; + const Type yi = *(this->library->target_values[i]); + + // Build design-row x̃_i = [1, lib[i][0], ..., lib[i][E-1]] + // We compute A and b without materializing x̃_i as a vector. + // x̃_i[0] = 1 (intercept), x̃_i[k] = lib[i][k-1] for k >= 1. + for (size_t k = 0; k < p; ++k) { + Type xk = (k == 0) ? Type(1) : *(this->library->embedded_values[i * E + k - 1]); + b_vec[k] += w * xk * yi; + for (size_t jj = 0; jj < p; ++jj) { + Type xj = (jj == 0) ? Type(1) : *(this->library->embedded_values[i * E + jj - 1]); + A[k * p + jj] += w * xk * xj; + } + } + } + + // ----------------------------------------------------------------------- + // Step 4: Solve A β = b via Cholesky decomposition (Eigen::LDLT). + // On exit b_vec holds the solution β. + // ----------------------------------------------------------------------- + Eigen::Matrix A_mat(p, p); + for (size_t i = 0; i < p; ++i) { + for (size_t j = 0; j < p; ++j) { + A_mat(i, j) = A[i * p + j]; + } + } + Eigen::Matrix b_mat(p); + for (size_t i = 0; i < p; ++i) { + b_mat(i) = b_vec[i]; + } + + Eigen::LDLT> ldlt(A_mat); + Eigen::Matrix sol = ldlt.solve(b_mat).eval(); + for (size_t i = 0; i < p; ++i) { + b_vec[i] = sol(i); + } + + // ----------------------------------------------------------------------- + // Step 5: Predict ŷ = [1, q[0], ..., q[E-1]] · β + // ----------------------------------------------------------------------- + Type prediction = b_vec[0]; // intercept + for (size_t j = 0; j < E; ++j) { + prediction += b_vec[j + 1] * query_point[j]; + } + return prediction; + } + + /** + * @brief Predict for all rows of a test embedding matrix. + * + * Iterates over every row of @p test_embedding, builds the query coordinate + * vector, calls predict_one(), and stores the results in `predictions`. + * After this call: predictions.size() == test_embedding.n_rows. + * + * @param test_embedding The query manifold. n_cols must equal + * embedding_dimension. + * @throws std::runtime_error / std::invalid_argument forwarded from + * predict_one() and ValidateLibrary(). + */ + virtual void predict( + const DelayEmbeddingMatrix& test_embedding) override { + this->ValidateLibrary(); + + const size_t E = this->embedding_dimension; + if (test_embedding.n_cols != E) { + throw std::invalid_argument( + "SMapProjection::predict: test embedding n_cols does not match " + "embedding_dimension."); + } + + const size_t n_test = test_embedding.n_rows; + this->predictions.resize(n_test); + + fims::Vector query_point(E); + for (size_t row = 0; row < n_test; ++row) { + for (size_t col = 0; col < E; ++col) { + query_point[col] = test_embedding.at(row, col); + } + this->predictions[row] = predict_one(query_point); + } + } +}; + +} // namespace fims_edm + +#endif // FIMS_EDM_SMAP_PROJECTION_HPP diff --git a/inst/include/edm/utilities/edm_distance_weights.hpp b/inst/include/edm/utilities/edm_distance_weights.hpp new file mode 100644 index 000000000..e8296b43b --- /dev/null +++ b/inst/include/edm/utilities/edm_distance_weights.hpp @@ -0,0 +1,259 @@ +/** + * @file edm_distance_weights.hpp + * @brief Modular distance and weighting utilities for EDM prediction algorithms. + * + * @details This header provides free functions in the fims_edm namespace that + * are shared across all EDM prediction algorithms: + * + * - SquaredEuclideanDistance(): pairwise squared distance between two + * embedding rows (raw-pointer interface for CppAD compatibility). + * - SimplexWeights(): exponential weights normalized by d_min (used by + * Simplex Projection, Sugihara & May 1990). + * - SMapWeights(): exponential or Gaussian weights scaled by theta. + * Kernel choice is controlled by the SMapKernel enum: + * - SMapKernel::kExponential (default): w = exp(-θ * d / d̄) + * Classic formulation from Sugihara (1994) and rEDM. + * - SMapKernel::kGaussian: w = exp(-θ² * (d/D)²) + * Gaussian kernel from Esguerra & Munch (2024); avoids sqrt inside + * the AD tape. + * - NormalizeWeights(): in-place normalization so weights sum to 1. + * + * ### Design rationale + * All weight vectors are passed by reference so callers can pre-allocate + * storage and avoid repeated heap allocation inside hot prediction loops. + * Raw-pointer distance computation keeps the inner loops visible to TMB's + * CppAD tracer, which cannot trace STL algorithms (e.g., std::transform, + * std::inner_product). + * + * @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_EDM_DISTANCE_WEIGHTS_HPP +#define FIMS_EDM_DISTANCE_WEIGHTS_HPP + +#include +#include +#include +#include + +#include "../../common/fims_math.hpp" + +namespace fims_edm { + +// --------------------------------------------------------------------------- +// Kernel selector +// --------------------------------------------------------------------------- + +/** + * @brief Selects the weighting kernel used by SMapWeights(). + * + * @details Both kernels are controlled by the nonlinearity parameter theta + * (θ ≥ 0). When θ = 0 every row has weight 1 regardless of the kernel + * chosen, so the distinction only matters for θ > 0. + * + * | Value | Formula | Source | + * |---|---|---| + * | kExponential | w = exp(−θ · d / d̄) | Sugihara (1994), rEDM | + * | kGaussian | w = exp(−θ² · (d/D)²) | Esguerra & Munch (2024) | + * + * The Gaussian kernel has slightly fatter tails and avoids a sqrt() call + * inside the AD tape (distances are kept squared throughout). + * Per Steve Munch (pers. comm.): "the difference in performance is rarely + * large; use whatever is easier/faster." + */ +enum class SMapKernel { + kExponential, ///< Classic exponential kernel: exp(-θ * d / d̄) + kGaussian ///< Gaussian kernel: exp(-θ² * (d/D)²) +}; + +// --------------------------------------------------------------------------- +// Distance +// --------------------------------------------------------------------------- + +/** + * @brief Compute the squared Euclidean distance between two embedding rows. + * + * @details Operates on raw pointer rows (as returned by the flat + * embedded_values storage in DelayEmbeddingMatrix) rather than on + * fims::Vector so the loop is visible to TMB's CppAD tracer. + * + * @tparam Type Numeric scalar type (double or TMB AD scalar). + * @param a_row Pointer to the first element of row a (length >= dim). + * @param b_row Pointer to the first element of row b (length >= dim). + * @param dim Number of coordinates per row (embedding dimension E). + * @return Squared Euclidean distance between the two rows. + * @throws std::invalid_argument if dim == 0. + */ +template +Type SquaredEuclideanDistance(const Type* a_row, const Type* b_row, + size_t dim) { + if (dim == 0) { + throw std::invalid_argument( + "SquaredEuclideanDistance: dim must be greater than 0."); + } + Type dist = Type(0); + for (size_t k = 0; k < dim; ++k) { + Type diff = a_row[k] - b_row[k]; + dist += diff * diff; + } + return dist; +} + +// --------------------------------------------------------------------------- +// Normalization +// --------------------------------------------------------------------------- + +/** + * @brief Normalize a weight vector in-place so that the weights sum to 1. + * + * @details All weights are divided by their total sum. If the sum is zero + * (e.g., all distances are infinite and all weights are zero), a + * std::runtime_error is thrown. + * + * @tparam Type Numeric scalar type. + * @param weights Weight vector to normalize in-place. + * @throws std::runtime_error if the weight sum is effectively zero. + */ +template +void NormalizeWeights(std::vector& weights) { + Type total = Type(0); + for (size_t i = 0; i < weights.size(); ++i) { + total += weights[i]; + } + if (total <= Type(0)) { + throw std::runtime_error( + "NormalizeWeights: weight sum is zero. " + "All distances may be infinite or weights all zero."); + } + for (size_t i = 0; i < weights.size(); ++i) { + weights[i] /= total; + } +} + +// --------------------------------------------------------------------------- +// Simplex Projection weights +// --------------------------------------------------------------------------- + +/** + * @brief Compute exponential weights for Simplex Projection. + * + * @details Implements the Sugihara & May (1990) exponential weighting scheme: + * \f[ + * w_i = \exp\!\Bigl(-\frac{d_i}{d_{\min} + \varepsilon}\Bigr) + * \f] + * where \f$d_{\min}\f$ is the smallest distance in @p sq_distances and + * \f$\varepsilon\f$ guards against division by zero when the query is an + * exact match. + * + * Weights are NOT normalized here; call NormalizeWeights() separately if + * required, or divide by their sum manually. + * + * @tparam Type Numeric scalar type. + * @param sq_distances Vector of squared Euclidean distances from the query + * point to each selected nearest neighbor. + * @param[out] weights Output weight vector. Resized to sq_distances.size(). + * @param epsilon Small constant added to d_min (default 1e-12). + * @throws std::invalid_argument if sq_distances is empty. + */ +template +void SimplexWeights(const std::vector& sq_distances, + std::vector& weights, double epsilon = 1e-12) { + if (sq_distances.empty()) { + throw std::invalid_argument("SimplexWeights: sq_distances must not be empty."); + } + + // Find d_min + Type d_min = sq_distances[0]; + for (size_t i = 1; i < sq_distances.size(); ++i) { + if (sq_distances[i] < d_min) { + d_min = sq_distances[i]; + } + } + + weights.resize(sq_distances.size()); + for (size_t i = 0; i < sq_distances.size(); ++i) { + Type exponent = sq_distances[i] / (d_min + Type(epsilon)); + weights[i] = fims_math::exp(-exponent); + } +} + +// --------------------------------------------------------------------------- +// S-Map weights +// --------------------------------------------------------------------------- + +/** + * @brief Compute S-Map weights with a user-selectable kernel. + * + * @details Two kernel choices are available via the SMapKernel enum: + * + * **Exponential** (SMapKernel::kExponential, default — Sugihara 1994 / rEDM): + * \f[ + * w_i = \exp\!\Bigl(-\theta \cdot \frac{d_i}{\bar{d}}\Bigr) + * \f] + * where \f$\bar{d}\f$ is the mean squared distance across ALL library rows. + * + * **Gaussian** (SMapKernel::kGaussian — Esguerra & Munch 2024): + * \f[ + * w_i = \exp\!\Bigl(-\theta^2 \cdot \Bigl(\frac{d_i}{D}\Bigr)^2\Bigr) + * \f] + * where \f$D\f$ is the mean squared distance (same denominator, kept squared + * to avoid a sqrt() inside the AD tape). + * + * For \f$\theta = 0\f$ both kernels reduce to \f$w_i = 1\f$ (global OLS). + * Per Munch (pers. comm.) the performance difference between kernels is + * rarely large; the Gaussian kernel avoids a sqrt() inside the AD tape. + * + * Weights are NOT normalized here; call NormalizeWeights() separately. + * + * @tparam Type Numeric scalar type (double or TMB AD scalar). + * @param sq_distances Squared Euclidean distances from the query point to + * every library row (length = n_library_rows). + * @param[out] weights Output weight vector, resized to sq_distances.size(). + * @param theta Nonlinearity parameter theta >= 0 (default 1.0). + * @param kernel Kernel choice: SMapKernel::kExponential (default) or + * SMapKernel::kGaussian. + * @param epsilon Small constant added to the scale denominator to guard + * against division by zero (default 1e-12). + * @throws std::invalid_argument if sq_distances is empty or theta < 0. + */ +template +void SMapWeights(const std::vector& sq_distances, + std::vector& weights, double theta = 1.0, + SMapKernel kernel = SMapKernel::kExponential, + double epsilon = 1e-12) { + if (sq_distances.empty()) { + throw std::invalid_argument("SMapWeights: sq_distances must not be empty."); + } + if (theta < 0.0) { + throw std::invalid_argument("SMapWeights: theta must be >= 0."); + } + + // Compute mean squared distance across all library rows (used by both + // kernels as the scale denominator D or d̄). + Type d_mean = Type(0); + for (size_t i = 0; i < sq_distances.size(); ++i) { + d_mean += sq_distances[i]; + } + d_mean /= Type(sq_distances.size()); + + weights.resize(sq_distances.size()); + for (size_t i = 0; i < sq_distances.size(); ++i) { + Type exponent; + if (kernel == SMapKernel::kGaussian) { + // Gaussian kernel: exp(-θ² · (d_i / D)²) + // d_i and D are already squared, so (d_i/D)² = d_i²/D². Keeping + // everything in squared-distance units avoids any sqrt() call. + Type ratio = sq_distances[i] / (d_mean + Type(epsilon)); + exponent = Type(theta * theta) * ratio * ratio; + } else { + // Exponential kernel (default): exp(-θ · d_i / d̄) + exponent = Type(theta) * sq_distances[i] / (d_mean + Type(epsilon)); + } + weights[i] = fims_math::exp(-exponent); + } +} + +} // namespace fims_edm + +#endif // FIMS_EDM_DISTANCE_WEIGHTS_HPP diff --git a/inst/include/edm/utilities/edm_kernels.hpp b/inst/include/edm/utilities/edm_kernels.hpp new file mode 100644 index 000000000..1dda4c9b7 --- /dev/null +++ b/inst/include/edm/utilities/edm_kernels.hpp @@ -0,0 +1,168 @@ +/** + * @file edm_kernels.hpp + * @brief Covariance kernel utilities for GP-EDM. + * + * @details Provides the Squared Exponential kernel with Automatic Relevance + * Determination (ARD) used by GPEdmProjection. The kernel follows the + * formulation of Munch et al. (2017) and the GPEDM R package by Tanya Rogers: + * + * \f[ + * K(x_i, x_j) = \sigma^2 \exp\!\Bigl( + * -\sum_{d=1}^{E} \phi_d (x_{id} - x_{jd})^2 + * \Bigr) + * \f] + * + * where: + * - \f$\phi_d \geq 0\f$ is the **inverse length-scale** for embedding + * dimension \f$d\f$ (ARD parameter): large \f$\phi_d\f$ means dimension + * \f$d\f$ matters a lot; \f$\phi_d \approx 0\f$ means it is irrelevant. + * - \f$\sigma^2 > 0\f$ is the signal variance (prior amplitude). + * + * The full covariance matrix used for prediction is: + * \f[ + * \Sigma = K + v_e I + * \f] + * where \f$v_e > 0\f$ is the process noise variance (nugget) added to the + * diagonal for numerical stability and to model observation uncertainty. + * + * All routines use raw pointer loops for CppAD traceability. + * + * @references + * Munch, S. B., Poynor, V., and Arriaza, J. L. (2017). Circumventing + * structural uncertainty: a Bayesian perspective on nonlinear forecasting + * for ecology. Ecological Complexity, 32: 134. + * + * Rogers, T. (2023). GPEDM: Gaussian process empirical dynamic modeling. + * https://github.com/tanyarogers/GPEDM + * + * @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_EDM_KERNELS_HPP +#define FIMS_EDM_KERNELS_HPP + +#include +#include +#include +#include + +#include "../../common/fims_math.hpp" + +namespace fims_edm { + +/** + * @brief Compute a single ARD squared-exponential kernel element. + * + * @details Evaluates + * \f[ + * k(x_i, x_j) = \sigma^2 \exp\!\Bigl(-\sum_{d=0}^{E-1} \phi_d + * (x_{id} - x_{jd})^2 \Bigr) + * \f] + * where \f$x_i\f$ and \f$x_j\f$ are embedding rows of length E, accessed via + * raw pointers. The summation is done in Type arithmetic so the result is + * differentiable under TMB/CppAD. + * + * @tparam Type Numeric scalar type (double or TMB AD scalar). + * @param row_i Pointer to the first element of embedding row i (length E). + * @param row_j Pointer to the first element of embedding row j (length E). + * @param E Embedding dimension (number of elements in each row). + * @param phi Per-dimension inverse length-scales (length E). Each entry + * must be >= 0. + * @param sigma2 Signal variance (must be > 0). + * @return \f$k(x_i, x_j)\f$ + */ +template +Type ARDKernelElement(const Type* row_i, const Type* row_j, size_t E, + const std::vector& phi, double sigma2) { + Type sq_dist = Type(0); + for (size_t d = 0; d < E; ++d) { + Type diff = row_i[d] - row_j[d]; + sq_dist += Type(phi[d]) * diff * diff; + } + return Type(sigma2) * fims_math::exp(-sq_dist); +} + +/** + * @brief Build the N×N ARD kernel matrix for the library. + * + * @details Fills the flat row-major vector @p K (length N*N) with: + * \f[ + * K[i \cdot N + j] = \sigma^2 \exp\!\Bigl(-\sum_d \phi_d + * (x_{id} - x_{jd})^2\Bigr) + * \f] + * and then adds the process noise @p ve to every diagonal element: + * \f[ + * K[i \cdot N + i] \mathrel{+}= v_e + * \f] + * + * The matrix is symmetric; only the upper triangle is computed and mirrored. + * + * @tparam Type Numeric scalar type (double or TMB AD scalar). + * @param embedded_values Flat row-major pointer array for the library + * (length N * E); element [row * E + col] points to + * the value at (row, col) of the embedding matrix. + * @param N Number of library rows. + * @param E Embedding dimension. + * @param phi Per-dimension inverse length-scales (length E). + * @param sigma2 Signal variance (> 0). + * @param ve Process noise variance added to the diagonal (> 0). + * @param[out] K Output flat row-major covariance matrix (N × N). + * Resized to N*N and overwritten. + */ +template +void BuildCovarianceMatrix( + const std::vector& embedded_values, size_t N, size_t E, + const std::vector& phi, double sigma2, double ve, + std::vector& K) { + K.resize(N * N); + + for (size_t i = 0; i < N; ++i) { + const Type* row_i = embedded_values[i * E]; + for (size_t j = i; j < N; ++j) { + const Type* row_j = embedded_values[j * E]; + Type k_ij = ARDKernelElement(row_i, row_j, E, phi, sigma2); + K[i * N + j] = k_ij; + K[j * N + i] = k_ij; // exploit symmetry + } + // Add nugget to diagonal + K[i * N + i] += Type(ve); + } +} + +/** + * @brief Compute the N-element covariance vector between a query point and + * every library row. + * + * @details Fills @p k_star so that: + * \f[ + * k^*[i] = \sigma^2 \exp\!\Bigl(-\sum_d \phi_d (q_d - x_{id})^2\Bigr) + * \f] + * + * Note: the nugget \f$v_e\f$ is NOT added here — it is only on the + * training covariance matrix diagonal. + * + * @tparam Type Numeric scalar type (double or TMB AD scalar). + * @param query Pointer to the E-element query row. + * @param embedded_values Flat row-major library pointer array (N * E). + * @param N Number of library rows. + * @param E Embedding dimension. + * @param phi Per-dimension inverse length-scales (length E). + * @param sigma2 Signal variance. + * @param[out] k_star Output vector of length N. Resized and overwritten. + */ +template +void BuildKStarVector(const Type* query, + const std::vector& embedded_values, + size_t N, size_t E, const std::vector& phi, + double sigma2, std::vector& k_star) { + k_star.resize(N); + for (size_t i = 0; i < N; ++i) { + const Type* row_i = embedded_values[i * E]; + k_star[i] = ARDKernelElement(query, row_i, E, phi, sigma2); + } +} + +} // namespace fims_edm + +#endif // FIMS_EDM_KERNELS_HPP diff --git a/inst/include/edm/utilities/edm_linear_algebra.hpp b/inst/include/edm/utilities/edm_linear_algebra.hpp new file mode 100644 index 000000000..00af420f2 --- /dev/null +++ b/inst/include/edm/utilities/edm_linear_algebra.hpp @@ -0,0 +1,23 @@ +/** + * @file edm_linear_algebra.hpp + * @brief Shared linear algebra utilities for EDM prediction algorithms. + * + * @details Linear solves across EDM functors (SMapProjection and GPEdmProjection) + * are standardized on Eigen::LDLT Cholesky factorization for numerical stability + * and differentiable linear algebra. + * + * @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_EDM_LINEAR_ALGEBRA_HPP +#define FIMS_EDM_LINEAR_ALGEBRA_HPP + +namespace fims_edm { + +// Linear algebra for EDM predictors (S-Map, GP-EDM) uses Eigen::LDLT directly +// for matrix factorization and solving linear systems. + +} // namespace fims_edm + +#endif // FIMS_EDM_LINEAR_ALGEBRA_HPP diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_edm.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_edm.hpp index 8e9a6eee6..9609716eb 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_edm.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_edm.hpp @@ -9,7 +9,7 @@ #ifndef FIMS_INTERFACE_RCPP_RCPP_OBJECTS_RCPP_EDM_HPP #define FIMS_INTERFACE_RCPP_RCPP_OBJECTS_RCPP_EDM_HPP -#include "../../../edm/functors/delay_embedding.hpp" +#include "../../../edm/edm.hpp" #include "rcpp_interface_base.hpp" /** @@ -427,4 +427,426 @@ class DelayEmbeddingInterface : public EDMInterfaceBase { #endif }; +/** + * @brief Helper to map RealVector memory from DelayEmbeddingInterface to a transient + * fims_edm::DelayEmbeddingMatrix representation for functors. + */ +inline fims_edm::DelayEmbeddingMatrix to_matrix(DelayEmbeddingInterface* de) { + fims_edm::DelayEmbeddingMatrix mat; + mat.n_rows = de->target_values.size(); + mat.n_cols = mat.n_rows > 0 ? (de->embedded_values.size() / mat.n_rows) : 0; + + if (de->embedded_values.size() > 0 && de->embedded_values.storage_m) { + mat.embedded_values.resize(de->embedded_values.size()); + for (size_t i = 0; i < de->embedded_values.size(); i++) { + mat.embedded_values[i] = &((*de->embedded_values.storage_m)[i]); + } + } + + if (de->target_values.size() > 0 && de->target_values.storage_m) { + mat.target_values.resize(de->target_values.size()); + for (size_t i = 0; i < de->target_values.size(); i++) { + mat.target_values[i] = &((*de->target_values.storage_m)[i]); + } + } + + if (de->embedded_uncertainty.size() > 0 && de->embedded_uncertainty.storage_m) { + mat.embedded_uncertainty.resize(de->embedded_uncertainty.size()); + for (size_t i = 0; i < de->embedded_uncertainty.size(); i++) { + mat.embedded_uncertainty[i] = &((*de->embedded_uncertainty.storage_m)[i]); + } + } + + if (de->target_uncertainty.size() > 0 && de->target_uncertainty.storage_m) { + mat.target_uncertainty.resize(de->target_uncertainty.size()); + for (size_t i = 0; i < de->target_uncertainty.size(); i++) { + mat.target_uncertainty[i] = &((*de->target_uncertainty.storage_m)[i]); + } + } + + return mat; +} + +/** + * @brief Rcpp interface for Simplex projection. + */ +class SimplexProjectionInterface : public EDMInterfaceBase { + public: + /** + * @brief The embedding dimension (E). + */ + uint32_t embedding_dimension; + /** + * @brief The number of neighbors (k). + */ + int32_t n_neighbors; + + /** + * @brief The constructor. + */ + SimplexProjectionInterface() : EDMInterfaceBase(), embedding_dimension(0), n_neighbors(0) { + EDMInterfaceBase::live_objects[this->id] = + std::make_shared(*this); + FIMSRcppInterfaceBase::fims_interface_objects.push_back( + EDMInterfaceBase::live_objects[this->id]); + } + + /** + * @brief Copy constructor. + */ + SimplexProjectionInterface(const SimplexProjectionInterface& other) + : EDMInterfaceBase(other), + embedding_dimension(other.embedding_dimension), + n_neighbors(other.n_neighbors) {} + + /** + * @brief The destructor. + */ + virtual ~SimplexProjectionInterface() {} + + /** + * @brief Get the ID of the object. + */ + virtual uint32_t get_id() { return this->id; } + + /** + * @brief Predict the values using Simplex projection. + * @param lib_id The ID of the library delay embedding. + * @param test_id The ID of the test delay embedding. + * @return A vector of predictions. + */ + Rcpp::NumericVector predict(uint32_t lib_id, uint32_t test_id) { + auto it_lib = EDMInterfaceBase::live_objects.find(lib_id); + auto it_test = EDMInterfaceBase::live_objects.find(test_id); + if (it_lib == EDMInterfaceBase::live_objects.end() || + it_test == EDMInterfaceBase::live_objects.end()) { + throw std::invalid_argument("SimplexProjectionInterface::predict: invalid lib_id or test_id."); + } + + DelayEmbeddingInterface* lib = dynamic_cast(it_lib->second.get()); + DelayEmbeddingInterface* test = dynamic_cast(it_test->second.get()); + if (!lib || !test) { + throw std::invalid_argument("SimplexProjectionInterface::predict: live objects are not of type DelayEmbeddingInterface."); + } + + fims_edm::DelayEmbeddingMatrix lib_matrix = to_matrix(lib); + fims_edm::DelayEmbeddingMatrix test_matrix = to_matrix(test); + + fims_edm::SimplexProjection sp; + sp.library = &lib_matrix; + sp.embedding_dimension = this->embedding_dimension; + sp.n_neighbors = this->n_neighbors; + + sp.predict(test_matrix); + + Rcpp::NumericVector ret(sp.predictions.size()); + for (size_t i = 0; i < sp.predictions.size(); ++i) { + ret[i] = sp.predictions[i]; + } + return ret; + } + + /** + * @brief Finalizes the interface object. + */ + virtual void finalize() { + this->finalized = true; + } + + /** + * @brief Converts to JSON representation. + */ + virtual std::string to_json() { + std::stringstream ss; + ss << "{\n"; + ss << " \"module_name\": \"EDM\",\n"; + ss << " \"module_type\": \"SimplexProjection\",\n"; + ss << " \"module_id\": " << this->id << ",\n"; + ss << " \"embedding_dimension\": " << this->embedding_dimension << ",\n"; + ss << " \"n_neighbors\": " << this->n_neighbors << "\n"; + ss << "}"; + return ss.str(); + } +}; + +/** + * @brief Rcpp interface for S-Map projection. + */ +class SMapProjectionInterface : public EDMInterfaceBase { + public: + /** + * @brief The embedding dimension (E). + */ + uint32_t embedding_dimension; + /** + * @brief The S-Map localization parameter (theta). + */ + double theta; + /** + * @brief The S-Map weighting kernel (exponential or gaussian). + */ + std::string kernel; + + /** + * @brief The constructor. + */ + SMapProjectionInterface() : EDMInterfaceBase(), embedding_dimension(0), theta(1.0), kernel("exponential") { + EDMInterfaceBase::live_objects[this->id] = + std::make_shared(*this); + FIMSRcppInterfaceBase::fims_interface_objects.push_back( + EDMInterfaceBase::live_objects[this->id]); + } + + /** + * @brief Copy constructor. + */ + SMapProjectionInterface(const SMapProjectionInterface& other) + : EDMInterfaceBase(other), + embedding_dimension(other.embedding_dimension), + theta(other.theta), + kernel(other.kernel) {} + + /** + * @brief The destructor. + */ + virtual ~SMapProjectionInterface() {} + + /** + * @brief Get the ID of the object. + */ + virtual uint32_t get_id() { return this->id; } + + /** + * @brief Predict the values using S-Map projection. + * @param lib_id The ID of the library delay embedding. + * @param test_id The ID of the test delay embedding. + * @return A vector of predictions. + */ + Rcpp::NumericVector predict(uint32_t lib_id, uint32_t test_id) { + auto it_lib = EDMInterfaceBase::live_objects.find(lib_id); + auto it_test = EDMInterfaceBase::live_objects.find(test_id); + if (it_lib == EDMInterfaceBase::live_objects.end() || + it_test == EDMInterfaceBase::live_objects.end()) { + throw std::invalid_argument("SMapProjectionInterface::predict: invalid lib_id or test_id."); + } + + DelayEmbeddingInterface* lib = dynamic_cast(it_lib->second.get()); + DelayEmbeddingInterface* test = dynamic_cast(it_test->second.get()); + if (!lib || !test) { + throw std::invalid_argument("SMapProjectionInterface::predict: live objects are not of type DelayEmbeddingInterface."); + } + + fims_edm::DelayEmbeddingMatrix lib_matrix = to_matrix(lib); + fims_edm::DelayEmbeddingMatrix test_matrix = to_matrix(test); + + fims_edm::SMapProjection sm; + sm.library = &lib_matrix; + sm.embedding_dimension = this->embedding_dimension; + sm.theta = this->theta; + if (this->kernel == "gaussian") { + sm.kernel = fims_edm::SMapKernel::kGaussian; + } else { + sm.kernel = fims_edm::SMapKernel::kExponential; + } + + sm.predict(test_matrix); + + Rcpp::NumericVector ret(sm.predictions.size()); + for (size_t i = 0; i < sm.predictions.size(); ++i) { + ret[i] = sm.predictions[i]; + } + return ret; + } + + /** + * @brief Finalizes the interface object. + */ + virtual void finalize() { + this->finalized = true; + } + + /** + * @brief Converts to JSON representation. + */ + virtual std::string to_json() { + std::stringstream ss; + ss << "{\n"; + ss << " \"module_name\": \"EDM\",\n"; + ss << " \"module_type\": \"SMapProjection\",\n"; + ss << " \"module_id\": " << this->id << ",\n"; + ss << " \"embedding_dimension\": " << this->embedding_dimension << ",\n"; + ss << " \"theta\": " << this->theta << ",\n"; + ss << " \"kernel\": \"" << this->kernel << "\"\n"; + ss << "}"; + return ss.str(); + } +}; + +/** + * @brief Rcpp interface for GP-EDM projection. + */ +class GPEdmProjectionInterface : public EDMInterfaceBase { + public: + /** + * @brief The embedding dimension (E). + */ + uint32_t embedding_dimension; + /** + * @brief The ARD length-scale parameters (phi). + */ + Rcpp::NumericVector phi; + /** + * @brief The signal variance parameter (sigma2). + */ + double sigma2; + /** + * @brief The observation noise variance/nugget parameter (ve). + */ + double ve; + + /** + * @brief The constructor. + */ + GPEdmProjectionInterface() : EDMInterfaceBase(), embedding_dimension(0), sigma2(1.0), ve(0.1) { + EDMInterfaceBase::live_objects[this->id] = + std::make_shared(*this); + FIMSRcppInterfaceBase::fims_interface_objects.push_back( + EDMInterfaceBase::live_objects[this->id]); + } + + /** + * @brief Copy constructor. + */ + GPEdmProjectionInterface(const GPEdmProjectionInterface& other) + : EDMInterfaceBase(other), + embedding_dimension(other.embedding_dimension), + phi(other.phi), + sigma2(other.sigma2), + ve(other.ve) {} + + /** + * @brief The destructor. + */ + virtual ~GPEdmProjectionInterface() {} + + /** + * @brief Get the ID of the object. + */ + virtual uint32_t get_id() { return this->id; } + + /** + * @brief Fits the optimized hyperparameters for GP-EDM. + * @param lib_id The ID of the library delay embedding. + * @return A list of optimized parameters. + */ + Rcpp::List fit(uint32_t lib_id) { + auto it_lib = EDMInterfaceBase::live_objects.find(lib_id); + if (it_lib == EDMInterfaceBase::live_objects.end()) { + throw std::invalid_argument("GPEdmProjectionInterface::fit: invalid lib_id."); + } + + DelayEmbeddingInterface* lib = dynamic_cast(it_lib->second.get()); + if (!lib) { + throw std::invalid_argument("GPEdmProjectionInterface::fit: live object is not of type DelayEmbeddingInterface."); + } + fims_edm::DelayEmbeddingMatrix lib_matrix = to_matrix(lib); + + fims_edm::GPEdmProjection gp; + gp.library = &lib_matrix; + gp.embedding_dimension = this->embedding_dimension; + gp.sigma2 = this->sigma2; + gp.ve = this->ve; + if (this->phi.size() > 0) { + gp.phi.resize(this->phi.size()); + for (int i = 0; i < this->phi.size(); ++i) { + gp.phi[i] = this->phi[i]; + } + } + + gp.fit(); + + this->sigma2 = gp.sigma2; + this->ve = gp.ve; + this->phi = Rcpp::NumericVector(gp.phi.size()); + for (size_t i = 0; i < gp.phi.size(); ++i) { + this->phi[i] = gp.phi[i]; + } + + return Rcpp::List::create( + Rcpp::Named("phi") = this->phi, + Rcpp::Named("sigma2") = this->sigma2, + Rcpp::Named("ve") = this->ve + ); + } + + /** + * @brief Predict the values using GP-EDM projection. + * @param lib_id The ID of the library delay embedding. + * @param test_id The ID of the test delay embedding. + * @return A vector of predictions. + */ + Rcpp::NumericVector predict(uint32_t lib_id, uint32_t test_id) { + auto it_lib = EDMInterfaceBase::live_objects.find(lib_id); + auto it_test = EDMInterfaceBase::live_objects.find(test_id); + if (it_lib == EDMInterfaceBase::live_objects.end() || + it_test == EDMInterfaceBase::live_objects.end()) { + throw std::invalid_argument("GPEdmProjectionInterface::predict: invalid lib_id or test_id."); + } + + DelayEmbeddingInterface* lib = dynamic_cast(it_lib->second.get()); + DelayEmbeddingInterface* test = dynamic_cast(it_test->second.get()); + if (!lib || !test) { + throw std::invalid_argument("GPEdmProjectionInterface::predict: live objects are not of type DelayEmbeddingInterface."); + } + + fims_edm::DelayEmbeddingMatrix lib_matrix = to_matrix(lib); + fims_edm::DelayEmbeddingMatrix test_matrix = to_matrix(test); + + fims_edm::GPEdmProjection gp; + gp.library = &lib_matrix; + gp.embedding_dimension = this->embedding_dimension; + gp.sigma2 = this->sigma2; + gp.ve = this->ve; + if (this->phi.size() > 0) { + gp.phi.resize(this->phi.size()); + for (int i = 0; i < this->phi.size(); ++i) { + gp.phi[i] = this->phi[i]; + } + } else { + gp.phi.assign(this->embedding_dimension, 0.1); + } + + gp.predict(test_matrix); + + Rcpp::NumericVector ret(gp.predictions.size()); + for (size_t i = 0; i < gp.predictions.size(); ++i) { + ret[i] = gp.predictions[i]; + } + return ret; + } + + /** + * @brief Finalizes the interface object. + */ + virtual void finalize() { + this->finalized = true; + } + + /** + * @brief Converts to JSON representation. + */ + virtual std::string to_json() { + std::stringstream ss; + ss << "{\n"; + ss << " \"module_name\": \"EDM\",\n"; + ss << " \"module_type\": \"GPEdmProjection\",\n"; + ss << " \"module_id\": " << this->id << ",\n"; + ss << " \"embedding_dimension\": " << this->embedding_dimension << ",\n"; + ss << " \"sigma2\": " << this->sigma2 << ",\n"; + ss << " \"ve\": " << this->ve << "\n"; + ss << "}"; + return ss.str(); + } +}; + #endif diff --git a/man/Cpp_classes.Rd b/man/Cpp_classes.Rd index 0f43cea8b..7fe04a4df 100644 --- a/man/Cpp_classes.Rd +++ b/man/Cpp_classes.Rd @@ -5,6 +5,7 @@ \alias{AgeComp} \alias{BevertonHoltRecruitment} \alias{CatchAtAge} +\alias{DelayEmbedding} \alias{DlnormDistribution} \alias{DmultinomDistribution} \alias{DnormDistribution} @@ -36,6 +37,7 @@ docs. \item \href{https://noaa-fims.github.io/FIMS/doxygen/classAgeCompDataInterface.html}{AgeComp} \item \href{https://noaa-fims.github.io/FIMS/doxygen/classBevertonHoltRecruitmentInterface.html}{BevertonHoltRecruitment} \item \href{https://noaa-fims.github.io/FIMS/doxygen/classCatchAtAgeInterface.html}{CatchAtAge} +\item \href{https://noaa-fims.github.io/FIMS/doxygen/classDelayEmbeddingInterface.html}{DelayEmbedding} \item \href{https://noaa-fims.github.io/FIMS/doxygen/classDlnormDistributionsInterface.html}{DlnormDistribution} \item \href{https://noaa-fims.github.io/FIMS/doxygen/classDmultinomDistributionsInterface.html}{DmultinomDistribution} \item \href{https://noaa-fims.github.io/FIMS/doxygen/classDnormDistributionsInterface.html}{DnormDistribution} diff --git a/man/FIMSFrame.Rd b/man/FIMSFrame.Rd index 6fbff87aa..3497ecc0a 100644 --- a/man/FIMSFrame.Rd +++ b/man/FIMSFrame.Rd @@ -16,7 +16,7 @@ See the \code{data_big} object in FIMS, e.g., \code{data(data_big, package = "FI An object of the S4 class \code{FIMSFrame} class, or one of its child classes, is validated and then returned. All objects will at a minimum have a slot called \code{data} to store the input data frame. Additional slots are dependent -on the child class. Use \code{\link[methods:showClass]{methods::showClass()}} to see all available slots. +on the child class. Use \code{\link[methods:RClassUtils]{methods::showClass()}} to see all available slots. } \description{ All constructor functions take a single input and build an object specific diff --git a/man/create_edm_embedding.Rd b/man/create_edm_embedding.Rd new file mode 100644 index 000000000..b12140083 --- /dev/null +++ b/man/create_edm_embedding.Rd @@ -0,0 +1,76 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/fimsframe.R +\name{create_edm_embedding} +\alias{create_edm_embedding} +\title{Create a delay embedding and store it inside a FIMSFrame object} +\usage{ +create_edm_embedding( + x, + series_type, + series_name, + E, + tau, + drop_missing = TRUE, + uncertainty_name = NULL, + embedding_name = paste0(series_type, "_", series_name, "_E", E, "_tau", tau) +) +} +\arguments{ +\item{x}{A \code{\link[=FIMSFrame]{FIMSFrame()}} object.} + +\item{series_type}{A single string giving the \code{type} value to filter on, +e.g. \code{"index"} or \code{"landings"}.} + +\item{series_name}{A single string giving the \code{name} (fleet / survey) value +to filter on, e.g. \code{"survey1"}.} + +\item{E}{A positive integer — the embedding dimension (number of lagged +coordinates per row).} + +\item{tau}{A positive integer — the time lag between successive coordinates.} + +\item{drop_missing}{Logical (default \code{TRUE}). When \code{TRUE}, rows containing +the FIMS missing-value sentinel \code{-999} are dropped from the embedding +matrix.} + +\item{uncertainty_name}{A single string giving the \code{name} column value for +an uncertainty series in the data slot (e.g. \code{"survey1_sd"}), or \code{NULL} +(default) to skip uncertainty propagation. The uncertainty series must +have the same length as the value series after filtering on +\code{series_type} and arranging by \code{timing}.} + +\item{embedding_name}{A single string used as the key under which the +result is stored in \code{edm_embeddings}. Defaults to +\code{paste0(series_type, "_", series_name, "_E", E, "_tau", tau)}.} +} +\value{ +A \code{\link[=FIMSFrame]{FIMSFrame()}} object identical to \code{x} with the new embedding +appended to the \code{edm_embeddings} slot. +} +\description{ +Pulls a univariate time series from the \code{data} slot of a \code{\link[=FIMSFrame]{FIMSFrame()}} +object, constructs the delay embedding matrix using the Rcpp +\code{DelayEmbedding} module, and returns a \strong{new} FIMSFrame with the result +stored in the \code{edm_embeddings} slot under \code{embedding_name}. +} +\details{ +The time series is extracted by filtering \code{get_data(x)} on \code{type} and +\code{name} (fleet / survey name). Values encoded as \code{-999} are treated as the +FIMS missing-data sentinel. When \code{drop_missing = TRUE}, +\code{construct_drop_missing()} is called so that rows containing \code{-999} are +removed from the embedding matrix; when \code{FALSE}, \code{construct()} is called +and those sentinel values remain in the matrix. + +Uncertainty propagation follows the same pointer-based pattern as the value +fields. When \code{uncertainty_name} is supplied, the function extracts a +matching time series (filtered on the same \code{series_type} and ordered by +\code{timing}) and passes it to the Rcpp module. The resulting embedding list +will include \code{embedded_uncertainty} and \code{target_uncertainty} numeric +vectors laid out identically to \code{embedded_values} and \code{target_values} +respectively. When \code{uncertainty_name} is \code{NULL} (the default) those fields +are absent from the stored list. + +The returned FIMSFrame is identical to \code{x} except that \code{edm_embeddings} +contains one additional element named \code{embedding_name}. +} +\keyword{FIMSFrame} diff --git a/man/data_big.Rd b/man/data_big.Rd index 2702a9621..c7d906249 100644 --- a/man/data_big.Rd +++ b/man/data_big.Rd @@ -9,7 +9,7 @@ A data frame with 2808 observations of 8 variables: \describe{ \item{type}{The type of input the row contains. Allowed types include -age_comp, age_to_length_conversion, index, landings, length_comp, and weight_at_age.} +age_comp, age_to_length_conversion, edm, index, landings, length_comp, and weight_at_age.} \item{fleet}{A character string providing the name of the information source that the input was collected from, e.g., \code{"Trawl fishery"} or \code{"age_zero_survey"}. Entries can be \code{NA} for information that is not diff --git a/man/get_FIMSFrame.Rd b/man/get_FIMSFrame.Rd index 025131f8a..c6643f8c1 100644 --- a/man/get_FIMSFrame.Rd +++ b/man/get_FIMSFrame.Rd @@ -32,6 +32,8 @@ \alias{get_n_lengths} \alias{get_n_lengths,FIMSFrame-method} \alias{get_n_lengths,data.frame-method} +\alias{get_edm_embeddings} +\alias{get_edm_embeddings,FIMSFrame-method} \title{Get a slot in a FIMSFrame object} \usage{ get_data(x) @@ -93,6 +95,10 @@ get_n_lengths(x) \S4method{get_n_lengths}{FIMSFrame}(x) \S4method{get_n_lengths}{data.frame}(x) + +get_edm_embeddings(x) + +\S4method{get_edm_embeddings}{FIMSFrame}(x) } \arguments{ \item{x}{An object returned from \code{\link[=FIMSFrame]{FIMSFrame()}}.} @@ -127,6 +133,14 @@ the model. \code{\link[=get_n_lengths]{get_n_lengths()}} returns an integer specifying the number of length bins used in the model. + +\code{\link[=get_edm_embeddings]{get_edm_embeddings()}} returns a named list of delay-embedding results +stored in the FIMSFrame object. Each element is itself a named list with +fields \code{series_name}, \code{E} (embedding dimension), \code{tau} (time lag), +\code{drop_missing} (logical), \code{n_rows}, \code{n_cols}, \code{embedded_values} (a numeric +vector of the flattened row-major matrix), and \code{target_values} (a numeric +vector of the target time-point values x_t for each row). +Returns an empty list when \code{\link[=create_edm_embedding]{create_edm_embedding()}} has not yet been called. } \description{ There is an accessor function for each slot in the S4 class \code{FIMSFrame}, diff --git a/man/lognormal.Rd b/man/lognormal.Rd index fbaa75048..1cc20277a 100644 --- a/man/lognormal.Rd +++ b/man/lognormal.Rd @@ -11,7 +11,7 @@ multinomial(link = "logit") } \arguments{ \item{link}{A string specifying the model link function. For example, -\code{"identity"} or \code{"log"} are appropriate names for the \code{\link[stats:gaussian]{stats::gaussian()}} +\code{"identity"} or \code{"log"} are appropriate names for the \code{\link[stats:family]{stats::gaussian()}} distribution. \code{"log"} and \code{"logit"} are the defaults for the lognormal and the multinomial, respectively.} } @@ -36,7 +36,7 @@ a_family[["link"]] \seealso{ \itemize{ \item \code{\link[stats:family]{stats::family()}} -\item \code{\link[stats:gaussian]{stats::gaussian()}} +\item \code{\link[stats:family]{stats::gaussian()}} \item \code{\link[stats:glm]{stats::glm()}} \item \code{\link[stats:power]{stats::power()}} \item \code{\link[stats:make.link]{stats::make.link()}} diff --git a/man/model_.Rd b/man/model_.Rd index 438442a02..6decbf5d9 100644 --- a/man/model_.Rd +++ b/man/model_.Rd @@ -20,6 +20,8 @@ \alias{model_age_to_length_conversion} \alias{model_age_to_length_conversion,FIMSFrame-method} \alias{model_age_to_length_conversion,data.frame-method} +\alias{model_edm_matrix} +\alias{model_edm_matrix,FIMSFrame-method} \title{Get a vector of data to be passed to a FIMS module from a FIMSFrame object} \usage{ model_landings(x, fleet) @@ -57,6 +59,10 @@ model_age_to_length_conversion(x) \S4method{model_age_to_length_conversion}{FIMSFrame}(x) \S4method{model_age_to_length_conversion}{data.frame}(x) + +model_edm_matrix(x, embedding_name) + +\S4method{model_edm_matrix}{FIMSFrame}(x, embedding_name) } \arguments{ \item{x}{An object returned from \code{\link[=FIMSFrame]{FIMSFrame()}}.} @@ -64,6 +70,9 @@ model_age_to_length_conversion(x) \item{fleet}{A string, or vector of strings, specifying the name of the fleet(s) of interest that you want data for. The strings must exactly match strings in the \code{"fleet"} column of \code{get_data(x)}.} + +\item{embedding_name}{A single string matching one of the names in +\code{\link[=get_edm_embeddings]{get_edm_embeddings()}}.} } \value{ All of the \verb{model_*()} functions return vectors of data. Currently, the @@ -125,5 +134,13 @@ multiple values are present across a single age and length, they are averaged because \code{age_to_length_conversion} data cannot vary across fleets or time. } + +\subsection{\code{model_edm_matrix()}}{ + +Returns a numeric matrix (n_rows × E) for a named delay embedding stored +inside the FIMSFrame object. The matrix can be passed directly to an EDM +prediction module. Rows correspond to embedded time points and columns +correspond to the lagged coordinates \link{x_t, x_{t-tau}, ..., x_{t-(E-1)tau}}. +} } \keyword{FIMSFrame} diff --git a/src/rcpp_edm.cpp b/src/rcpp_edm.cpp index b595c3322..1fa893ca6 100644 --- a/src/rcpp_edm.cpp +++ b/src/rcpp_edm.cpp @@ -43,4 +43,35 @@ void register_edm(Rcpp::Module& m) { .method("construct_drop_missing_with_uncertainty", &DelayEmbeddingInterface::construct_drop_missing_with_uncertainty) .method("at", &DelayEmbeddingInterface::at); + + Rcpp::class_( + "SimplexProjection", + "Exposes Simplex projection functor to R.") + .constructor() + .field("embedding_dimension", &SimplexProjectionInterface::embedding_dimension) + .field("n_neighbors", &SimplexProjectionInterface::n_neighbors) + .method("get_id", &SimplexProjectionInterface::get_id) + .method("predict", &SimplexProjectionInterface::predict); + + Rcpp::class_( + "SMapProjection", + "Exposes S-Map projection functor to R.") + .constructor() + .field("embedding_dimension", &SMapProjectionInterface::embedding_dimension) + .field("theta", &SMapProjectionInterface::theta) + .field("kernel", &SMapProjectionInterface::kernel) + .method("get_id", &SMapProjectionInterface::get_id) + .method("predict", &SMapProjectionInterface::predict); + + Rcpp::class_( + "GPEdmProjection", + "Exposes GP-EDM projection functor to R.") + .constructor() + .field("embedding_dimension", &GPEdmProjectionInterface::embedding_dimension) + .field("phi", &GPEdmProjectionInterface::phi) + .field("sigma2", &GPEdmProjectionInterface::sigma2) + .field("ve", &GPEdmProjectionInterface::ve) + .method("get_id", &GPEdmProjectionInterface::get_id) + .method("fit", &GPEdmProjectionInterface::fit) + .method("predict", &GPEdmProjectionInterface::predict); } diff --git a/tests/gtest/CMakeLists.txt b/tests/gtest/CMakeLists.txt index 2fe3c52cc..25bc23115 100644 --- a/tests/gtest/CMakeLists.txt +++ b/tests/gtest/CMakeLists.txt @@ -283,6 +283,42 @@ target_link_libraries(edm_delay_embedding gtest_discover_tests(edm_delay_embedding) +# test_edm_simplex_projection.cpp +add_executable(edm_simplex_projection + test_edm_simplex_projection.cpp +) + +target_link_libraries(edm_simplex_projection + gtest_main + fims_test +) + +gtest_discover_tests(edm_simplex_projection) + +# test_edm_smap_projection.cpp +add_executable(edm_smap_projection + test_edm_smap_projection.cpp +) + +target_link_libraries(edm_smap_projection + gtest_main + fims_test +) + +gtest_discover_tests(edm_smap_projection) + +# test_edm_gp_projection.cpp +add_executable(edm_gp_projection + test_edm_gp_projection.cpp +) + +target_link_libraries(edm_gp_projection + gtest_main + fims_test +) + +gtest_discover_tests(edm_gp_projection) + # test_population_CatchAtAge_CalculateUnfishedInitial.cpp add_executable(population_CatchAtAge_CalculateUnfishedInitial test_population_CatchAtAge_CalculateUnfishedInitial.cpp @@ -348,3 +384,26 @@ target_link_libraries(def_FIMSLog_clear fims_test ) gtest_discover_tests(def_FIMSLog_clear) + +# test_distributions_mvnorm.cpp +add_executable(distributions_mvnorm + test_distributions_mvnorm.cpp +) +add_as_invoker_manifest(distributions_mvnorm) +target_link_libraries(distributions_mvnorm + gtest_main + fims_test +) +gtest_discover_tests(distributions_mvnorm) + +# test_distributions_beta.cpp +add_executable(distributions_beta + test_distributions_beta.cpp +) +add_as_invoker_manifest(distributions_beta) +target_link_libraries(distributions_beta + gtest_main + fims_test +) +gtest_discover_tests(distributions_beta) + diff --git a/tests/gtest/test_distributions_beta.cpp b/tests/gtest/test_distributions_beta.cpp new file mode 100644 index 000000000..5ec93bfb4 --- /dev/null +++ b/tests/gtest/test_distributions_beta.cpp @@ -0,0 +1,74 @@ +/** + * @file test_distributions_beta.cpp + * @brief GoogleTest unit tests for BetaLPDF distribution functor. + */ + +#include + +#include +#include + +#include "distributions/functors/beta_lpdf.hpp" + +TEST(BetaLPDF, StandardUniformBetaLogDensity) { + // Beta(1, 1) on (0, 1) is Uniform(0, 1) -> pdf = 1, log_pdf = 0. + fims_distributions::BetaLPDF beta; + beta.shape1.resize(1); beta.shape1[0] = 1.0; + beta.shape2.resize(1); beta.shape2[0] = 1.0; + + beta.observed_values.resize(1); + beta.observed_values[0] = 0.5; + + beta.expected_values.resize(1); + beta.expected_values[0] = 0.0; + + beta.input_type = "prior"; + beta.priors.clear(); + beta.priors.push_back(&(beta.observed_values)); + + double lpdf = beta.evaluate(); + EXPECT_NEAR(lpdf, 0.0, 1e-6); +} + +TEST(BetaLPDF, SymmetricBetaLogDensity) { + // Beta(2, 2) on (0, 1): pdf(x) = 6 * x * (1 - x). At x = 0.5 -> pdf = 1.5 -> log_pdf = log(1.5) + fims_distributions::BetaLPDF beta; + beta.shape1.resize(1); beta.shape1[0] = 2.0; + beta.shape2.resize(1); beta.shape2[0] = 2.0; + + beta.observed_values.resize(1); + beta.observed_values[0] = 0.5; + + beta.expected_values.resize(1); + beta.expected_values[0] = 0.0; + + beta.input_type = "prior"; + beta.priors.clear(); + beta.priors.push_back(&(beta.observed_values)); + + double lpdf = beta.evaluate(); + double expected_lpdf = std::log(1.5); + EXPECT_NEAR(lpdf, expected_lpdf, 1e-6); +} + +TEST(BetaLPDF, ScaledBetaLogDensity) { + // Beta(2, 2) on (0, max=10): at x = 5.0 -> u = 0.5 -> pdf(x) = 1.5 / 10 = 0.15 -> log_pdf = log(0.15) + fims_distributions::BetaLPDF beta; + beta.shape1.resize(1); beta.shape1[0] = 2.0; + beta.shape2.resize(1); beta.shape2[0] = 2.0; + beta.scale.resize(1); beta.scale[0] = 10.0; + + beta.observed_values.resize(1); + beta.observed_values[0] = 5.0; + + beta.expected_values.resize(1); + beta.expected_values[0] = 0.0; + + beta.input_type = "prior"; + beta.priors.clear(); + beta.priors.push_back(&(beta.observed_values)); + + double lpdf = beta.evaluate(); + double expected_lpdf = std::log(0.15); + EXPECT_NEAR(lpdf, expected_lpdf, 1e-6); +} diff --git a/tests/gtest/test_distributions_mvnorm.cpp b/tests/gtest/test_distributions_mvnorm.cpp new file mode 100644 index 000000000..7ec7f37ae --- /dev/null +++ b/tests/gtest/test_distributions_mvnorm.cpp @@ -0,0 +1,71 @@ +/** + * @file test_distributions_mvnorm.cpp + * @brief GoogleTest unit tests for MVNormLPDF distribution functor. + */ + +#include + +#include +#include + +#include "distributions/functors/mvnorm_lpdf.hpp" + +TEST(MVNormLPDF, BivariateStandardNormalLogDensity) { + // Test bivariate standard normal distribution with mean (0, 0) and identity covariance + fims_distributions::MVNormLPDF mvnorm; + mvnorm.k_dim = 2; + + mvnorm.sigma_mat.resize(4); + mvnorm.sigma_mat[0] = 1.0; mvnorm.sigma_mat[1] = 0.0; + mvnorm.sigma_mat[2] = 0.0; mvnorm.sigma_mat[3] = 1.0; + + mvnorm.observed_values.resize(2); + mvnorm.observed_values[0] = 0.0; + mvnorm.observed_values[1] = 0.0; + + mvnorm.expected_values.resize(2); + mvnorm.expected_values[0] = 0.0; + mvnorm.expected_values[1] = 0.0; + + mvnorm.input_type = "prior"; + mvnorm.priors.clear(); + mvnorm.priors.push_back(&(mvnorm.observed_values)); + + double lpdf = mvnorm.evaluate(); + + // Analytical bivariate standard normal at (0, 0): + // -k/2 * log(2*pi) - 0.5 * log(1) - 0.5 * 0 = -log(2*pi) ≈ -1.837877 + double expected_lpdf = -std::log(2.0 * M_PI); + EXPECT_NEAR(lpdf, expected_lpdf, 1e-6); +} + +TEST(MVNormLPDF, CorrelatedBivariateNormalLogDensity) { + // Test bivariate normal with covariance Sigma = [[2, 1], [1, 2]] + // det(Sigma) = 3 + fims_distributions::MVNormLPDF mvnorm; + mvnorm.k_dim = 2; + + mvnorm.sigma_mat.resize(4); + mvnorm.sigma_mat[0] = 2.0; mvnorm.sigma_mat[1] = 1.0; + mvnorm.sigma_mat[2] = 1.0; mvnorm.sigma_mat[3] = 2.0; + + mvnorm.observed_values.resize(2); + mvnorm.observed_values[0] = 1.0; + mvnorm.observed_values[1] = 1.0; + + mvnorm.expected_values.resize(2); + mvnorm.expected_values[0] = 0.0; + mvnorm.expected_values[1] = 0.0; + + mvnorm.input_type = "prior"; + mvnorm.priors.clear(); + mvnorm.priors.push_back(&(mvnorm.observed_values)); + + double lpdf = mvnorm.evaluate(); + + // Quad form: x = [1, 1], Sigma^-1 = 1/3 * [[2, -1], [-1, 2]] + // x^T Sigma^-1 x = 1/3 * (2 - 1 - 1 + 2) = 2/3 + // Expected lpdf = -0.5 * (2/3 + log(3) + 2*log(2*pi)) + double expected_lpdf = -0.5 * (2.0 / 3.0 + std::log(3.0) + 2.0 * std::log(2.0 * M_PI)); + EXPECT_NEAR(lpdf, expected_lpdf, 1e-5); +} diff --git a/tests/gtest/test_edm_delay_embedding.cpp b/tests/gtest/test_edm_delay_embedding.cpp index 38e3ed786..409a93d87 100644 --- a/tests/gtest/test_edm_delay_embedding.cpp +++ b/tests/gtest/test_edm_delay_embedding.cpp @@ -1,58 +1,78 @@ #include "edm/functors/delay_embedding.hpp" #include "gtest/gtest.h" +// ### Layout note (h=1 default) +// With forecast_horizon=1 (default): +// n_rows = series.size() - lag_span - 1 +// embedded[row] = [series[t], series[t-tau], ...] where t = row + lag_span +// target[row] = series[t + 1] (one step ahead) + namespace { TEST(DelayEmbedding, BuildsEmbeddingWithUnitLag) { + // series={1,2,3,4,5}, E=3, tau=1, h=1 + // lag_span = 2, n_rows = 5 - 2 - 1 = 2 + // row 0: t=2 → embedded=[3,2,1], target=series[3]=4 + // row 1: t=3 → embedded=[4,3,2], target=series[4]=5 fims::Vector series = {1.0, 2.0, 3.0, 4.0, 5.0}; fims_edm::DelayEmbeddingMatrix embedding = fims_edm::MakeDelayEmbedding(series, 3, 1); - EXPECT_EQ(embedding.n_rows, 3); + EXPECT_EQ(embedding.n_rows, 2); EXPECT_EQ(embedding.n_cols, 3); - // target_values[row] points to x_t for that row - EXPECT_DOUBLE_EQ(*embedding.target_values[0], 3.0); // row 0 → x_2 = 3.0 - EXPECT_DOUBLE_EQ(*embedding.target_values[2], 5.0); // row 2 → x_4 = 5.0 + // target_values[row] = x_{t+1} + EXPECT_DOUBLE_EQ(*embedding.target_values[0], 4.0); // row 0 target = x_3 + EXPECT_DOUBLE_EQ(*embedding.target_values[1], 5.0); // row 1 target = x_4 // embedded_values: [x_t, x_{t-1}, x_{t-2}] EXPECT_DOUBLE_EQ(embedding.at(0, 0), 3.0); EXPECT_DOUBLE_EQ(embedding.at(0, 1), 2.0); EXPECT_DOUBLE_EQ(embedding.at(0, 2), 1.0); - EXPECT_DOUBLE_EQ(embedding.at(2, 0), 5.0); - EXPECT_DOUBLE_EQ(embedding.at(2, 1), 4.0); - EXPECT_DOUBLE_EQ(embedding.at(2, 2), 3.0); + EXPECT_DOUBLE_EQ(embedding.at(1, 0), 4.0); + EXPECT_DOUBLE_EQ(embedding.at(1, 1), 3.0); + EXPECT_DOUBLE_EQ(embedding.at(1, 2), 2.0); } TEST(DelayEmbedding, BuildsEmbeddingWithLargerLag) { + // series={10,20,30,40,50,60}, E=2, tau=2, h=1 + // lag_span = 2, n_rows = 6 - 2 - 1 = 3 + // row 0: t=2 → embedded=[30,10], target=series[3]=40 + // row 2: t=4 → embedded=[50,30], target=series[5]=60 fims::Vector series = {10.0, 20.0, 30.0, 40.0, 50.0, 60.0}; fims_edm::DelayEmbeddingMatrix embedding = fims_edm::MakeDelayEmbedding(series, 2, 2); - EXPECT_EQ(embedding.n_rows, 4); + EXPECT_EQ(embedding.n_rows, 3); EXPECT_EQ(embedding.n_cols, 2); - // Row 0: target is x_2 = 30.0 (first valid index given lag_span = 2) - EXPECT_DOUBLE_EQ(*embedding.target_values[0], 30.0); + // Row 0: t=2, target = x_3 = 40 + EXPECT_DOUBLE_EQ(*embedding.target_values[0], 40.0); EXPECT_DOUBLE_EQ(embedding.at(0, 0), 30.0); EXPECT_DOUBLE_EQ(embedding.at(0, 1), 10.0); - EXPECT_DOUBLE_EQ(embedding.at(3, 0), 60.0); - EXPECT_DOUBLE_EQ(embedding.at(3, 1), 40.0); + EXPECT_DOUBLE_EQ(embedding.at(2, 0), 50.0); + EXPECT_DOUBLE_EQ(embedding.at(2, 1), 30.0); } TEST(DelayEmbedding, SupportsEmbeddingDimensionOne) { + // series={4,5,6}, E=1, tau=1, h=1 + // lag_span=0, n_rows = 3 - 0 - 1 = 2 + // row 0: t=0, embedded=[4], target=5 + // row 1: t=1, embedded=[5], target=6 fims::Vector series = {4.0, 5.0, 6.0}; fims_edm::DelayEmbeddingMatrix embedding = fims_edm::MakeDelayEmbedding(series, 1, 1); - EXPECT_EQ(embedding.n_rows, 3); + EXPECT_EQ(embedding.n_rows, 2); EXPECT_EQ(embedding.n_cols, 1); EXPECT_DOUBLE_EQ(embedding.at(0, 0), 4.0); - EXPECT_DOUBLE_EQ(embedding.at(2, 0), 6.0); + EXPECT_DOUBLE_EQ(embedding.at(1, 0), 5.0); + EXPECT_DOUBLE_EQ(*embedding.target_values[0], 5.0); + EXPECT_DOUBLE_EQ(*embedding.target_values[1], 6.0); } TEST(DelayEmbedding, RejectsInvalidInputs) { @@ -62,30 +82,34 @@ TEST(DelayEmbedding, RejectsInvalidInputs) { std::invalid_argument); EXPECT_THROW(fims_edm::MakeDelayEmbedding(series, 2, 0), std::invalid_argument); + // E=4, tau=1: lag_span=3, need size >= 3+1+1=5, only have 3 EXPECT_THROW(fims_edm::MakeDelayEmbedding(series, 4, 1), std::invalid_argument); + // E=2, tau=3: lag_span=3, need size >= 3+1+1=5, only have 3 EXPECT_THROW(fims_edm::MakeDelayEmbedding(series, 2, 3), std::invalid_argument); } TEST(DelayEmbedding, DropsRowsWithMissingValues) { + // series={1,-999,3,4,5}, E=2, tau=1, h=1 + // Full embedding (before drop): n_rows = 5-1-1 = 3 + // row 0: t=1, embedded=[-999,1], target=3 → has missing → drop + // row 1: t=2, embedded=[3,-999], target=4 → has missing → drop + // row 2: t=3, embedded=[4,3], target=5 → valid + // After drop: 1 valid row const double missing_value = -999.0; - fims::Vector series = {1.0, 2.0, missing_value, 4.0, 5.0}; + fims::Vector series = {1.0, missing_value, 3.0, 4.0, 5.0}; fims_edm::DelayEmbeddingMatrix embedding = fims_edm::MakeDelayEmbeddingDropMissing(series, 2, 1, missing_value); - EXPECT_EQ(embedding.n_rows, 2); + EXPECT_EQ(embedding.n_rows, 1); EXPECT_EQ(embedding.n_cols, 2); - // Row 0: target was x_1 = 2.0; Row 1: target was x_4 = 5.0 - EXPECT_DOUBLE_EQ(*embedding.target_values[0], 2.0); - EXPECT_DOUBLE_EQ(*embedding.target_values[1], 5.0); - - EXPECT_DOUBLE_EQ(embedding.at(0, 0), 2.0); - EXPECT_DOUBLE_EQ(embedding.at(0, 1), 1.0); - EXPECT_DOUBLE_EQ(embedding.at(1, 0), 5.0); - EXPECT_DOUBLE_EQ(embedding.at(1, 1), 4.0); + // Retained row: embedded=[4,3], target=5 + EXPECT_DOUBLE_EQ(*embedding.target_values[0], 5.0); + EXPECT_DOUBLE_EQ(embedding.at(0, 0), 4.0); + EXPECT_DOUBLE_EQ(embedding.at(0, 1), 3.0); } TEST(DelayEmbedding, ThrowsOnOutOfBoundsMatrixAccess) { @@ -98,54 +122,45 @@ TEST(DelayEmbedding, ThrowsOnOutOfBoundsMatrixAccess) { } TEST(DelayEmbedding, DynamicallyReflectsOriginalSeriesModifications) { + // series={1,2,3,4,5}, E=3, tau=1, h=1 → n_rows=2 + // row 0: t=2, embedded=[3,2,1], target=4 fims::Vector series = {1.0, 2.0, 3.0, 4.0, 5.0}; fims_edm::DelayEmbeddingMatrix embedding = fims_edm::MakeDelayEmbedding(series, 3, 1); - // Initial values EXPECT_DOUBLE_EQ(embedding.at(0, 0), 3.0); - EXPECT_DOUBLE_EQ(embedding.at(2, 0), 5.0); - EXPECT_DOUBLE_EQ(*embedding.target_values[0], 3.0); + EXPECT_DOUBLE_EQ(*embedding.target_values[0], 4.0); - // Modify the original series elements - series[2] = 99.0; - series[4] = 100.0; + // Modify the original series → pointers reflect change immediately + series[3] = 99.0; // target of row 0 - // Both embedded_values and target_values dynamically reflect changes - EXPECT_DOUBLE_EQ(embedding.at(0, 0), 99.0); - EXPECT_DOUBLE_EQ(embedding.at(2, 0), 100.0); EXPECT_DOUBLE_EQ(*embedding.target_values[0], 99.0); - EXPECT_DOUBLE_EQ(*embedding.target_values[2], 100.0); } TEST(DelayEmbedding, PropagatesUncertaintyVectors) { - fims::Vector series = {1.0, 2.0, 3.0, 4.0, 5.0}; - // Uncertainty series: sigma_i = 0.1 * (i + 1) + // series={1,2,3,4,5}, E=3, tau=1, h=1 → n_rows=2 + // uncertainty={0.1,0.2,0.3,0.4,0.5} + // row 0: t=2, embedded uncertainty=[0.3,0.2,0.1], target uncertainty=0.4 + // row 1: t=3, embedded uncertainty=[0.4,0.3,0.2], target uncertainty=0.5 + fims::Vector series = {1.0, 2.0, 3.0, 4.0, 5.0}; fims::Vector uncertainty = {0.1, 0.2, 0.3, 0.4, 0.5}; fims_edm::DelayEmbeddingMatrix embedding = fims_edm::MakeDelayEmbedding(series, 3, 1, uncertainty); - // Same shape as value embedding - EXPECT_EQ(embedding.n_rows, 3); + EXPECT_EQ(embedding.n_rows, 2); EXPECT_EQ(embedding.n_cols, 3); - // target_uncertainty[row] mirrors target_values[row] in the uncertainty space - // Row 0: target index = 2 → sigma_2 = 0.3 - EXPECT_DOUBLE_EQ(*embedding.target_uncertainty[0], 0.3); - // Row 2: target index = 4 → sigma_4 = 0.5 - EXPECT_DOUBLE_EQ(*embedding.target_uncertainty[2], 0.5); + // target_uncertainty[row] = sigma_{t+1} + EXPECT_DOUBLE_EQ(*embedding.target_uncertainty[0], 0.4); // sigma_3 + EXPECT_DOUBLE_EQ(*embedding.target_uncertainty[1], 0.5); // sigma_4 - // embedded_uncertainty layout mirrors embedded_values: - // Row 0, col 0: sigma_{t} = sigma_2 = 0.3 - // Row 0, col 1: sigma_{t-1} = sigma_1 = 0.2 - // Row 0, col 2: sigma_{t-2} = sigma_0 = 0.1 + // embedded_uncertainty layout mirrors embedded_values EXPECT_DOUBLE_EQ(*embedding.embedded_uncertainty[0 * 3 + 0], 0.3); EXPECT_DOUBLE_EQ(*embedding.embedded_uncertainty[0 * 3 + 1], 0.2); EXPECT_DOUBLE_EQ(*embedding.embedded_uncertainty[0 * 3 + 2], 0.1); - // Row 2, col 0: sigma_{t} = sigma_4 = 0.5 - EXPECT_DOUBLE_EQ(*embedding.embedded_uncertainty[2 * 3 + 0], 0.5); + EXPECT_DOUBLE_EQ(*embedding.embedded_uncertainty[1 * 3 + 0], 0.4); } TEST(DelayEmbedding, UncertaintyFieldsEmptyWhenNotProvided) { @@ -154,48 +169,45 @@ TEST(DelayEmbedding, UncertaintyFieldsEmptyWhenNotProvided) { fims_edm::DelayEmbeddingMatrix embedding = fims_edm::MakeDelayEmbedding(series, 3, 1); - // No uncertainty series supplied → both fields must be empty EXPECT_TRUE(embedding.embedded_uncertainty.empty()); EXPECT_TRUE(embedding.target_uncertainty.empty()); } TEST(DelayEmbedding, UncertaintyDynamicallyReflectsSeriesModifications) { - fims::Vector series = {1.0, 2.0, 3.0, 4.0, 5.0}; + // row 0: t=2, target uncertainty = sigma_{3} = 0.4 + fims::Vector series = {1.0, 2.0, 3.0, 4.0, 5.0}; fims::Vector uncertainty = {0.1, 0.2, 0.3, 0.4, 0.5}; fims_edm::DelayEmbeddingMatrix embedding = fims_edm::MakeDelayEmbedding(series, 3, 1, uncertainty); - // Initial uncertainty at target row 0 = 0.3 - EXPECT_DOUBLE_EQ(*embedding.target_uncertainty[0], 0.3); + EXPECT_DOUBLE_EQ(*embedding.target_uncertainty[0], 0.4); // Modify the original uncertainty series element - uncertainty[2] = 9.9; + uncertainty[3] = 9.9; - // Pointer-based: embedding reflects the update without copying EXPECT_DOUBLE_EQ(*embedding.target_uncertainty[0], 9.9); - EXPECT_DOUBLE_EQ(*embedding.embedded_uncertainty[0 * 3 + 0], 9.9); } TEST(DelayEmbedding, DropMissingPropagatesUncertainty) { + // series={1,-999,3,4,5}, E=2, tau=1, h=1 + // Only row 2 (t=3, embedded=[4,3], target=5) is valid after drop const double missing_value = -999.0; - fims::Vector series = {1.0, 2.0, missing_value, 4.0, 5.0}; - fims::Vector uncertainty = {0.1, 0.2, 0.3, 0.4, 0.5}; + fims::Vector series = {1.0, missing_value, 3.0, 4.0, 5.0}; + fims::Vector uncertainty = {0.1, 0.2, 0.3, 0.4, 0.5}; fims_edm::DelayEmbeddingMatrix embedding = fims_edm::MakeDelayEmbeddingDropMissing(series, 2, 1, missing_value, uncertainty); - // Two valid rows: row 0 → target x_1 = 2.0, row 1 → target x_4 = 5.0 - EXPECT_EQ(embedding.n_rows, 2); + EXPECT_EQ(embedding.n_rows, 1); - // target_uncertainty matches the retained rows - EXPECT_DOUBLE_EQ(*embedding.target_uncertainty[0], 0.2); // sigma_1 - EXPECT_DOUBLE_EQ(*embedding.target_uncertainty[1], 0.5); // sigma_4 + // target_uncertainty = sigma_{t+1} = sigma_4 = 0.5 + EXPECT_DOUBLE_EQ(*embedding.target_uncertainty[0], 0.5); - // embedded_uncertainty at row 0: [sigma_1, sigma_0] = [0.2, 0.1] - EXPECT_DOUBLE_EQ(*embedding.embedded_uncertainty[0 * 2 + 0], 0.2); - EXPECT_DOUBLE_EQ(*embedding.embedded_uncertainty[0 * 2 + 1], 0.1); + // embedded_uncertainty at row 0: [sigma_3, sigma_2] = [0.4, 0.3] + EXPECT_DOUBLE_EQ(*embedding.embedded_uncertainty[0 * 2 + 0], 0.4); + EXPECT_DOUBLE_EQ(*embedding.embedded_uncertainty[0 * 2 + 1], 0.3); } } // namespace diff --git a/tests/gtest/test_edm_gp_projection.cpp b/tests/gtest/test_edm_gp_projection.cpp new file mode 100644 index 000000000..7880662a7 --- /dev/null +++ b/tests/gtest/test_edm_gp_projection.cpp @@ -0,0 +1,412 @@ +/** + * @file test_edm_gp_projection.cpp + * @brief GoogleTest unit tests for GPEdmProjection. + * + * Tests cover: + * - Guard conditions (null library, wrong dimension, too-small library) + * - ARD kernel properties (flat kernel when phi=0; localized when phi large) + * - Prediction output size + * - Near-interpolation with near-zero nugget + * - Prediction changes with forecast horizon + * - fit() reduces the negative log posterior + * - Logistic map prediction stays in plausible range + * + * @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. + */ +#include +#include + +#include "gtest/gtest.h" +#include "../../inst/include/edm/edm.hpp" + +// --------------------------------------------------------------------------- +// Helpers +// --------------------------------------------------------------------------- + +/// Build a fims::Vector from a std::initializer_list. +static fims::Vector make_vec(std::initializer_list vals) { + fims::Vector v; + v.resize(vals.size()); + size_t i = 0; + for (double x : vals) v[i++] = x; + return v; +} + +/// Build a simple linear series 0, 1, 2, ..., n-1. +static fims::Vector linear_series(size_t n) { + fims::Vector s; + s.resize(n); + for (size_t i = 0; i < n; ++i) s[i] = static_cast(i); + return s; +} + +/// Logistic map: x_{t+1} = r * x_t * (1 - x_t), r=3.8, x0=0.4. +static fims::Vector logistic_map(size_t n, double r = 3.8, + double x0 = 0.4) { + fims::Vector s; + s.resize(n); + s[0] = x0; + for (size_t i = 1; i < n; ++i) s[i] = r * s[i - 1] * (1.0 - s[i - 1]); + return s; +} + +/// Set up a GPEdmProjection with library from the first n_lib rows of the +/// given embedding, E=E, phi all equal to phi_val. +struct GPSetup { + fims::Vector series; + fims_edm::DelayEmbeddingMatrix lib_emb; + fims_edm::DelayEmbeddingMatrix test_emb; + fims_edm::GPEdmProjection gp; + + GPSetup(fims::Vector s, size_t E, double phi_val = 0.5, + double sigma2_val = 1.0, double ve_val = 1e-4, + size_t forecast_horizon = 1) + : series(std::move(s)) { + lib_emb = fims_edm::MakeDelayEmbedding(series, E, /*tau=*/1, + fims::Vector(), + forecast_horizon); + test_emb = lib_emb; // use same embedding as test for predict() + + gp.library = &lib_emb; + gp.embedding_dimension = E; + gp.phi.assign(E, phi_val); + gp.sigma2 = sigma2_val; + gp.ve = ve_val; + } +}; + +// --------------------------------------------------------------------------- +// Guard / error tests +// --------------------------------------------------------------------------- + +TEST(GPEdmProjection, ThrowsOnNullLibrary) { + fims_edm::GPEdmProjection gp; + gp.embedding_dimension = 2; + gp.phi = {0.5, 0.5}; + + fims::Vector q = make_vec({1.0, 2.0}); + EXPECT_THROW(gp.predict_one(q), std::runtime_error); +} + +TEST(GPEdmProjection, ThrowsOnNullLibraryFit) { + fims_edm::GPEdmProjection gp; + gp.embedding_dimension = 2; + EXPECT_THROW(gp.fit(), std::runtime_error); +} + +TEST(GPEdmProjection, ThrowsWhenLibraryTooSmall) { + // Need at least 2 rows; E=2 needs series length >= 4 + // Use a 3-element series with E=2, horizon=1 → 1 row only + fims::Vector s = make_vec({1.0, 2.0, 3.0}); + fims_edm::DelayEmbeddingMatrix emb = + fims_edm::MakeDelayEmbedding(s, /*E=*/2, /*tau=*/1); + // emb has 1 row → too small + + fims_edm::GPEdmProjection gp; + gp.library = &emb; + gp.embedding_dimension = 2; + gp.phi = {0.5, 0.5}; + + fims::Vector q = make_vec({1.0, 2.0}); + EXPECT_THROW(gp.predict_one(q), std::runtime_error); + EXPECT_THROW(gp.fit(), std::runtime_error); +} + +TEST(GPEdmProjection, ThrowsOnQueryDimensionMismatch) { + auto s = linear_series(15); + GPSetup setup(s, /*E=*/2); + + // query has 3 elements, embedding_dimension = 2 + fims::Vector q = make_vec({1.0, 2.0, 3.0}); + EXPECT_THROW(setup.gp.predict_one(q), std::invalid_argument); +} + +TEST(GPEdmProjection, ThrowsOnPhiSizeMismatch) { + auto s = linear_series(15); + GPSetup setup(s, /*E=*/2); + setup.gp.phi = {0.5}; // wrong: should be length 2 + + fims::Vector q = make_vec({1.0, 2.0}); + EXPECT_THROW(setup.gp.predict_one(q), std::invalid_argument); +} + +TEST(GPEdmProjection, ThrowsOnTestEmbeddingDimensionMismatch) { + auto s = linear_series(15); + GPSetup setup(s, /*E=*/2); + + // Build a test embedding with E=3 + fims_edm::DelayEmbeddingMatrix bad_test = + fims_edm::MakeDelayEmbedding(s, /*E=*/3, /*tau=*/1); + EXPECT_THROW(setup.gp.predict(bad_test), std::invalid_argument); +} + +// --------------------------------------------------------------------------- +// Prediction output size +// --------------------------------------------------------------------------- + +TEST(GPEdmProjection, PredictFillsCorrectOutputSize) { + auto s = linear_series(20); + GPSetup setup(s, /*E=*/2); + + setup.gp.predict(setup.test_emb); + EXPECT_EQ(setup.gp.predictions.size(), setup.test_emb.n_rows); +} + +TEST(GPEdmProjection, PredictOneReturnsScalar) { + auto s = linear_series(20); + GPSetup setup(s, /*E=*/2); + + fims::Vector q = make_vec({5.0, 4.0}); + double pred = setup.gp.predict_one(q); + // Just check it's a finite number + EXPECT_TRUE(std::isfinite(pred)); +} + +// --------------------------------------------------------------------------- +// ARD kernel property tests +// --------------------------------------------------------------------------- + +TEST(GPEdmProjection, ZeroPhiGivesConstantKernel) { + // When all phi[d] = 0, exp(-0) = 1, so K[i,j] = sigma2 for all i != j. + // All library points contribute equally regardless of distance. + // Prediction must lie between min and max of targets. + auto s = logistic_map(30); + GPSetup setup(s, /*E=*/2, /*phi_val=*/0.0, /*sigma2=*/1.0, /*ve=*/1e-3); + + fims::Vector q = make_vec({0.5, 0.5}); + double pred = setup.gp.predict_one(q); + EXPECT_TRUE(std::isfinite(pred)); +} + +TEST(GPEdmProjection, LargePhiLocalizesKernel) { + // Direct verification of the ARD kernel property: + // K(x_i, x_j) = sigma2 * exp(-sum_d phi_d * (x_id - x_jd)^2) + // Larger phi => kernel decays faster with distance => smaller K for same gap. + const size_t E = 2; + const double sigma2 = 1.0; + + // Two library rows: row_a (near query), row_b (far from query) + std::vector row_a_vals = {0.5, 0.4}; + std::vector row_b_vals = {0.5, 0.5}; // small gap from row_a + std::vector ptrs = {row_a_vals.data(), row_b_vals.data()}; + + // Dummy Type pointers for ARDKernelElement (uses raw Type* not vector) + // We call the utility function directly. + std::vector phi_small(E, 0.1); + std::vector phi_large(E, 10.0); + + double k_small = fims_edm::ARDKernelElement( + row_a_vals.data(), row_b_vals.data(), E, phi_small, sigma2); + double k_large = fims_edm::ARDKernelElement( + row_a_vals.data(), row_b_vals.data(), E, phi_large, sigma2); + + // K(a, b) with large phi must be strictly smaller (faster decay). + // rows differ only on dim 1: gap = 0.5 - 0.4 = 0.1, sq = 0.01 + // K_small = exp(-0.1 * 0.01) = exp(-0.001) + // K_large = exp(-10.0 * 0.01) = exp(-0.1) + EXPECT_GT(k_small, k_large) + << "Larger phi should produce smaller cross-covariance (faster decay)."; + EXPECT_NEAR(k_small, std::exp(-0.1 * 0.01), 1e-10); + EXPECT_NEAR(k_large, std::exp(-10.0 * 0.01), 1e-10); +} + +TEST(GPEdmProjection, LargePhiGivesDifferentPredThanSmallPhi) { + // ARD sensitivity: changing phi should change prediction. + auto s = logistic_map(30); + const size_t E = 2; + fims_edm::DelayEmbeddingMatrix emb = + fims_edm::MakeDelayEmbedding(s, E, /*tau=*/1); + + fims::Vector q = make_vec({0.6, 0.3}); + + fims_edm::GPEdmProjection gp_small; + gp_small.library = &emb; + gp_small.embedding_dimension = E; + gp_small.phi = {0.1, 0.1}; + gp_small.sigma2 = 1.0; + gp_small.ve = 1e-3; + + fims_edm::GPEdmProjection gp_large; + gp_large.library = &emb; + gp_large.embedding_dimension = E; + gp_large.phi = {10.0, 10.0}; + gp_large.sigma2 = 1.0; + gp_large.ve = 1e-3; + + double pred_small = gp_small.predict_one(q); + double pred_large = gp_large.predict_one(q); + + EXPECT_NE(pred_small, pred_large); + EXPECT_TRUE(std::isfinite(pred_small)); + EXPECT_TRUE(std::isfinite(pred_large)); +} + +// --------------------------------------------------------------------------- +// Forecast horizon test +// --------------------------------------------------------------------------- + +TEST(GPEdmProjection, ForecastHorizonTwoGivesDifferentPredictions) { + auto s = logistic_map(40); + const size_t E = 2; + + fims_edm::DelayEmbeddingMatrix emb_h1 = + fims_edm::MakeDelayEmbedding(s, E, /*tau=*/1, fims::Vector(), + /*forecast_horizon=*/1); + fims_edm::DelayEmbeddingMatrix emb_h2 = + fims_edm::MakeDelayEmbedding(s, E, /*tau=*/1, fims::Vector(), + /*forecast_horizon=*/2); + + fims_edm::GPEdmProjection gp_h1, gp_h2; + gp_h1.library = &emb_h1; gp_h1.embedding_dimension = E; + gp_h1.phi = {0.5, 0.5}; gp_h1.sigma2 = 1.0; gp_h1.ve = 1e-3; + + gp_h2.library = &emb_h2; gp_h2.embedding_dimension = E; + gp_h2.phi = {0.5, 0.5}; gp_h2.sigma2 = 1.0; gp_h2.ve = 1e-3; + + fims::Vector q = make_vec({0.5, 0.4}); + double pred_h1 = gp_h1.predict_one(q); + double pred_h2 = gp_h2.predict_one(q); + + EXPECT_NE(pred_h1, pred_h2); + EXPECT_TRUE(std::isfinite(pred_h1)); + EXPECT_TRUE(std::isfinite(pred_h2)); +} + +// --------------------------------------------------------------------------- +// fit() tests +// --------------------------------------------------------------------------- + +TEST(GPEdmProjection, FitProducesFiniteHyperparameters) { + auto s = logistic_map(40); + const size_t E = 2; + fims_edm::DelayEmbeddingMatrix emb = + fims_edm::MakeDelayEmbedding(s, E, /*tau=*/1); + + fims_edm::GPEdmProjection gp; + gp.library = &emb; + gp.embedding_dimension = E; + // phi not set — fit() should initialize to 0.1 + + EXPECT_NO_THROW(gp.fit()); + + ASSERT_EQ(gp.phi.size(), E); + for (size_t d = 0; d < E; ++d) { + EXPECT_TRUE(std::isfinite(gp.phi[d])); + EXPECT_GE(gp.phi[d], 0.0); + } + EXPECT_TRUE(std::isfinite(gp.sigma2)); + EXPECT_GT(gp.sigma2, 0.0); + EXPECT_TRUE(std::isfinite(gp.ve)); + EXPECT_GT(gp.ve, 0.0); +} + +TEST(GPEdmProjection, FitFollowedByPredictGivesFiniteResults) { + auto s = logistic_map(50); + const size_t E = 2; + fims_edm::DelayEmbeddingMatrix lib_emb = + fims_edm::MakeDelayEmbedding(s, E, /*tau=*/1); + + fims_edm::GPEdmProjection gp; + gp.library = &lib_emb; + gp.embedding_dimension = E; + + gp.fit(); + gp.predict(lib_emb); + + ASSERT_EQ(gp.predictions.size(), lib_emb.n_rows); + for (size_t i = 0; i < gp.predictions.size(); ++i) { + EXPECT_TRUE(std::isfinite(gp.predictions[i])) + << "prediction[" << i << "] is not finite"; + } +} + +TEST(GPEdmProjection, FitInitializesPhiWhenEmpty) { + auto s = logistic_map(30); + const size_t E = 3; + fims_edm::DelayEmbeddingMatrix emb = + fims_edm::MakeDelayEmbedding(s, E, /*tau=*/1); + + fims_edm::GPEdmProjection gp; + gp.library = &emb; + gp.embedding_dimension = E; + // phi intentionally left empty + + gp.fit(); + EXPECT_EQ(gp.phi.size(), E); +} + +// --------------------------------------------------------------------------- +// Logistic map sanity check +// --------------------------------------------------------------------------- + +TEST(GPEdmProjection, LogisticMapPredictionInRange) { + // After fitting on logistic map data, predictions should be near [0,1] + auto s = logistic_map(60); + const size_t E = 2; + fims_edm::DelayEmbeddingMatrix lib_emb = + fims_edm::MakeDelayEmbedding(s, E, /*tau=*/1); + + fims_edm::GPEdmProjection gp; + gp.library = &lib_emb; + gp.embedding_dimension = E; + gp.fit(); + gp.predict(lib_emb); + + // In-sample predictions on a logistic map should be within [-0.5, 1.5] + // (GP posterior mean can slightly exceed [0,1] but should be close) + for (size_t i = 0; i < gp.predictions.size(); ++i) { + EXPECT_GT(gp.predictions[i], -0.5) + << "prediction[" << i << "] = " << gp.predictions[i]; + EXPECT_LT(gp.predictions[i], 1.5) + << "prediction[" << i << "] = " << gp.predictions[i]; + } +} + +// --------------------------------------------------------------------------- +// SMap still works after GaussianElimination was moved +// --------------------------------------------------------------------------- + +TEST(GPEdmProjection, SMapStillWorksAfterRefactor) { + // Regression test: SMapProjection now uses Eigen::LDLT for Cholesky decomposition. + // Verify it still produces correct output. + auto s = logistic_map(30); + const size_t E = 2; + fims_edm::DelayEmbeddingMatrix emb = + fims_edm::MakeDelayEmbedding(s, E, /*tau=*/1); + + fims_edm::SMapProjection smap; + smap.library = &emb; + smap.embedding_dimension = E; + smap.theta = 1.0; + + EXPECT_NO_THROW(smap.predict(emb)); + ASSERT_EQ(smap.predictions.size(), emb.n_rows); + for (size_t i = 0; i < smap.predictions.size(); ++i) { + EXPECT_TRUE(std::isfinite(smap.predictions[i])); + } +} + +TEST(GPEdmProjection, IllConditionedNearZeroNuggetHandledByLDLT) { + // Steve Munch / Andrea review: verify near-zero nugget (ve ~ 0) handles + // potential PSD / near-singular matrix using Eigen::LDLT without throwing or NaN. + auto s = logistic_map(30); + const size_t E = 2; + fims_edm::DelayEmbeddingMatrix emb = + fims_edm::MakeDelayEmbedding(s, E, /*tau=*/1); + + fims_edm::GPEdmProjection gp; + gp.library = &emb; + gp.embedding_dimension = E; + gp.phi = std::vector(E, 100.0); // large phi + gp.ve = 1e-6; // near-zero nugget + + EXPECT_NO_THROW({ + gp.predict(emb); + for (size_t i = 0; i < gp.predictions.size(); ++i) { + EXPECT_TRUE(std::isfinite(gp.predictions[i])); + } + }); +} + diff --git a/tests/gtest/test_edm_simplex_projection.cpp b/tests/gtest/test_edm_simplex_projection.cpp new file mode 100644 index 000000000..22d8733ca --- /dev/null +++ b/tests/gtest/test_edm_simplex_projection.cpp @@ -0,0 +1,225 @@ +#include "edm/edm.hpp" +#include "gtest/gtest.h" + +namespace { + +// --------------------------------------------------------------------------- +// Test 1: Default n_neighbors is 0 (predicts using E+1 internally) +// --------------------------------------------------------------------------- +TEST(SimplexProjection, DefaultNeighborCountIsZero) { + fims_edm::SimplexProjection sp; + sp.embedding_dimension = 3; + sp.time_lag = 1; + + // n_neighbors == 0 is the sentinel meaning "use E+1 internally". + EXPECT_EQ(sp.n_neighbors, 0u); + + // A library with 5 rows (≥ E+1 = 4) must succeed with n_neighbors = 0. + fims::Vector series = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0}; + fims_edm::DelayEmbeddingMatrix lib = + fims_edm::MakeDelayEmbedding(series, 3, 1); + + sp.library = &lib; + fims::Vector q = {5.0, 4.0, 3.0}; + EXPECT_NO_THROW(sp.predict_one(q)); +} + +// --------------------------------------------------------------------------- +// Test 2: Exact-match point dominates exponential weighting +// +// For the linear series x_t = t (E=1, tau=1, k=2, h=1): +// - Library rows: embedded={t}, target=x_{t+1}=t+1 +// series={1,2,3,4,5,6} → rows: [1]→2, [2]→3, [3]→4, [4]→5 +// - Query = {4.0} → exact match to embedded=4, target=5 +// - d_min = 0 → exponent for exact match = 0 → w=1 +// - exponent for second neighbor = large → w≈0 +// - Prediction ≈ target of exact match = 5.0 +// --------------------------------------------------------------------------- +TEST(SimplexProjection, ExactMatchDominatesWeighting) { + fims::Vector series = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + fims_edm::DelayEmbeddingMatrix lib = + fims_edm::MakeDelayEmbedding(series, 1, 1); // h=1 default + + fims_edm::SimplexProjection sp; + sp.embedding_dimension = 1; + sp.time_lag = 1; + sp.n_neighbors = 2; + sp.library = &lib; + + // Query matches embedded=4 exactly. Its target is 5 (one step ahead). + fims::Vector q = {4.0}; + double pred = sp.predict_one(q); + + // The exact-match row dominates → prediction ≈ 5.0 + EXPECT_NEAR(pred, 5.0, 1e-6); +} + +// --------------------------------------------------------------------------- +// Test 3: Uniform weights when all neighbors are equidistant +// +// Manually construct a library with 4 rows all at squared distance 1.0 +// from the query point (1.0, 0.0): +// Row 0: embedded=(2.0, 0.0), target=10.0 → d²=1 +// Row 1: embedded=(0.0, 0.0), target=20.0 → d²=1 +// Row 2: embedded=(1.0, 1.0), target=30.0 → d²=1 +// Row 3: embedded=(1.0,-1.0), target=40.0 → d²=1 +// All exponents = 1/(1+eps) ≈ 1 → equal weights → mean = 25.0. +// --------------------------------------------------------------------------- +TEST(SimplexProjection, AllEqualDistancesGiveUniformWeights) { + fims_edm::DelayEmbeddingMatrix lib; + lib.n_rows = 4; + lib.n_cols = 2; + + // Static storage so pointers remain valid for the test lifetime. + static double vals[8] = {2.0, 0.0, // row 0 + 0.0, 0.0, // row 1 + 1.0, 1.0, // row 2 + 1.0, -1.0}; // row 3 + static double tgts[4] = {10.0, 20.0, 30.0, 40.0}; + + lib.embedded_values.resize(8); + lib.target_values.resize(4); + for (size_t i = 0; i < 8; ++i) lib.embedded_values[i] = &vals[i]; + for (size_t i = 0; i < 4; ++i) lib.target_values[i] = &tgts[i]; + + fims_edm::SimplexProjection sp; + sp.embedding_dimension = 2; + sp.time_lag = 1; + sp.n_neighbors = 4; + sp.library = &lib; + + fims::Vector q = {1.0, 0.0}; + double pred = sp.predict_one(q); + + // All weights identical → prediction = arithmetic mean of targets. + EXPECT_NEAR(pred, 25.0, 1e-6); +} + +// --------------------------------------------------------------------------- +// Test 4: predict() fills predictions with n_test elements +// --------------------------------------------------------------------------- +TEST(SimplexProjection, PredictFillsOutputVector) { + fims::Vector series = {1.0, 2.0, 3.0, 4.0, 5.0, + 6.0, 7.0, 8.0, 9.0, 10.0}; + + fims_edm::DelayEmbeddingMatrix lib = + fims_edm::MakeDelayEmbedding(series, 2, 1); // h=1: n_rows = 10-1-1 = 8 + ASSERT_EQ(lib.n_rows, 8u); + + fims_edm::DelayEmbeddingMatrix test_emb = + fims_edm::MakeDelayEmbedding(series, 2, 1); + + fims_edm::SimplexProjection sp; + sp.embedding_dimension = 2; + sp.time_lag = 1; + sp.n_neighbors = 3; + sp.library = &lib; + + sp.predict(test_emb); + + EXPECT_EQ(sp.predictions.size(), test_emb.n_rows); +} + +// --------------------------------------------------------------------------- +// Test 5: Throws when library pointer is null +// --------------------------------------------------------------------------- +TEST(SimplexProjection, ThrowsOnNullLibrary) { + fims_edm::SimplexProjection sp; + sp.embedding_dimension = 2; + sp.library = nullptr; + + fims::Vector q = {1.0, 2.0}; + EXPECT_THROW(sp.predict_one(q), std::runtime_error); +} + +// --------------------------------------------------------------------------- +// Test 6: Throws when query dimension mismatches embedding_dimension +// --------------------------------------------------------------------------- +TEST(SimplexProjection, ThrowsOnQueryDimensionMismatch) { + fims::Vector series = {1.0, 2.0, 3.0, 4.0, 5.0}; + fims_edm::DelayEmbeddingMatrix lib = + fims_edm::MakeDelayEmbedding(series, 2, 1); + + fims_edm::SimplexProjection sp; + sp.embedding_dimension = 2; + sp.library = &lib; + + // Query has 3 elements but embedding_dimension is 2. + fims::Vector q_wrong = {1.0, 2.0, 3.0}; + EXPECT_THROW(sp.predict_one(q_wrong), std::invalid_argument); +} + +// --------------------------------------------------------------------------- +// Test 7: Throws when library has fewer rows than requested k +// --------------------------------------------------------------------------- +TEST(SimplexProjection, ThrowsWhenLibraryTooSmall) { + fims::Vector series = {1.0, 2.0, 3.0, 4.0}; + // E=2, tau=1, h=1 → n_rows = 4 - 1 - 1 = 2. + fims_edm::DelayEmbeddingMatrix lib = + fims_edm::MakeDelayEmbedding(series, 2, 1); + ASSERT_EQ(lib.n_rows, 2u); + + fims_edm::SimplexProjection sp; + sp.embedding_dimension = 2; + sp.library = &lib; + sp.n_neighbors = 10; // More neighbors than library rows. + + fims::Vector q = {2.0, 1.0}; + EXPECT_THROW(sp.predict_one(q), std::runtime_error); +} + +// --------------------------------------------------------------------------- +// Test 8: Logistic-map prediction stays within [0, 1] +// +// The logistic map x_{t+1} = r * x_t * (1 - x_t) with r=3.8 is bounded +// in [0, 1]. A reasonable Simplex prediction on the embedded series should +// also fall within that range. +// --------------------------------------------------------------------------- +TEST(SimplexProjection, LogisticMapPredictionInRange) { + const size_t N = 50; + const double r = 3.8; + fims::Vector series; + series.emplace_back(0.4); + for (size_t i = 1; i < N; ++i) { + double prev = series[i - 1]; + series.emplace_back(r * prev * (1.0 - prev)); + } + + fims_edm::DelayEmbeddingMatrix lib = + fims_edm::MakeDelayEmbedding(series, 3, 1); + + fims_edm::SimplexProjection sp; + sp.embedding_dimension = 3; + sp.time_lag = 1; + // n_neighbors = 0 → uses E+1 = 4 internally. + sp.library = &lib; + + // Query: the last three values of the series. + fims::Vector q = {series[N - 1], series[N - 2], series[N - 3]}; + double pred = sp.predict_one(q); + + // Logistic map is bounded in (0, 1). + EXPECT_GE(pred, 0.0); + EXPECT_LE(pred, 1.0); +} + +// --------------------------------------------------------------------------- +// Test 9: predict() throws if test embedding has wrong dimension +// --------------------------------------------------------------------------- +TEST(SimplexProjection, PredictThrowsOnTestDimensionMismatch) { + fims::Vector series = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + fims_edm::DelayEmbeddingMatrix lib = + fims_edm::MakeDelayEmbedding(series, 2, 1); + + // Test embedding uses E=3, but predictor expects E=2. + fims_edm::DelayEmbeddingMatrix test_emb = + fims_edm::MakeDelayEmbedding(series, 3, 1); + + fims_edm::SimplexProjection sp; + sp.embedding_dimension = 2; + sp.library = &lib; + + EXPECT_THROW(sp.predict(test_emb), std::invalid_argument); +} + +} // namespace diff --git a/tests/gtest/test_edm_smap_projection.cpp b/tests/gtest/test_edm_smap_projection.cpp new file mode 100644 index 000000000..889f75e48 --- /dev/null +++ b/tests/gtest/test_edm_smap_projection.cpp @@ -0,0 +1,326 @@ +// Instructions ---- +// This file follows the format used for EDM gtest files in FIMS. +// Tests cover: IO correctness, edge handling, and error handling. +// +// ### Data-layout note (post mentor feedback, 2024-06) +// MakeDelayEmbedding now takes a forecast_horizon parameter (default = 1). +// With the default: +// +// embedded_values[row][col] = &series[target_index - col * tau] +// target_values[row] = &series[target_index + 1] (x_{t+1}) +// +// This matches the standard EDM one-step-ahead formulation (Sugihara 1994; +// Esguerra & Munch 2024, Section 2.1): predict x_{t+1} from x_t. +// +// Consequence: the design matrix and target are now DIFFERENT. For +// arithmetic series {a, a+k, a+2k, ...} the OLS still fits a perfect line, +// but the slope now captures the actual +k step, not a degenerate self-fit. +// Arithmetic series still produce a singular design matrix for E>=2 because +// col-1 = col-0 - k (linear combo of intercept and col-0). All tests with +// E>=2 therefore use quadratic or logistic-map series. +// +// ### Kernel note (post mentor feedback, 2024-06) +// SMapProjection now exposes a `kernel` field: +// SMapKernel::kExponential (default): w = exp(-theta * d / d_mean) [rEDM] +// SMapKernel::kGaussian: w = exp(-theta^2 * (d/D)^2) [E&M 2024] + +#include "edm/edm.hpp" +#include "gtest/gtest.h" + +using namespace fims_edm; + +namespace { + +// --------------------------------------------------------------------------- +// Helper: build a DelayEmbeddingMatrix from a flat time series. +// Uses forecast_horizon = 1 (the new default) throughout. +// --------------------------------------------------------------------------- +DelayEmbeddingMatrix MakeLib(const fims::Vector& series, + size_t E, size_t tau = 1, + size_t h = 1) { + return MakeDelayEmbedding(series, E, tau, fims::Vector(), h); +} + +// --------------------------------------------------------------------------- +// Helper: generate the logistic map x_{t+1} = r * x_t * (1 - x_t) +// --------------------------------------------------------------------------- +fims::Vector LogisticMap(double x0, double r, size_t n) { + fims::Vector s(n); + s[0] = x0; + for (size_t i = 1; i < n; ++i) { + s[i] = r * s[i - 1] * (1.0 - s[i - 1]); + } + return s; +} + +// =========================================================================== +// Test 1: theta=0, E=1, arithmetic series → target = x_{t+1} = x_t + 1. +// +// Series: {1, 2, 3, ..., 10}. +// With h=1 and E=1: embedding = [x_t], target = x_{t+1} = x_t + 1. +// OLS fit: y = 1 + 1*x_t. Query q={7} → prediction ≈ 8. +// =========================================================================== +TEST(SMapProjection, ThetaZeroE1LinearSeriesOneStepAhead) { + fims::Vector series = {1.0, 2.0, 3.0, 4.0, 5.0, + 6.0, 7.0, 8.0, 9.0, 10.0}; + auto lib = MakeLib(series, 1); // h=1 default + + SMapProjection smap; + smap.library = &lib; + smap.embedding_dimension = 1; + smap.theta = 0.0; + + // Global OLS: y = 1 + x_t → prediction for x_t=7 is 8 + fims::Vector q = {7.0}; + double pred = smap.predict_one(q); + EXPECT_NEAR(pred, 8.0, 0.5); +} + +// =========================================================================== +// Test 2: theta=0, quadratic series with E=2. +// +// Series: {1, 4, 9, 16, 25, 36, 49, 64, 81}. +// With h=1, E=2: embedding = [x_t, x_{t-1}], target = x_{t+1}. +// The step pattern is non-constant so the design matrix is full-rank. +// For theta=0 (global OLS) the prediction should be in a plausible range. +// =========================================================================== +TEST(SMapProjection, ThetaZeroQuadraticE2FullRankDesignMatrix) { + fims::Vector series = {1.0, 4.0, 9.0, 16.0, 25.0, + 36.0, 49.0, 64.0, 81.0, 100.0}; + auto lib = MakeLib(series, 2); + + SMapProjection smap; + smap.library = &lib; + smap.embedding_dimension = 2; + smap.theta = 0.0; + + // Query [64, 49]: the global OLS should predict a value close to 81 + fims::Vector q = {64.0, 49.0}; + double pred = smap.predict_one(q); + EXPECT_NEAR(pred, 81.0, 5.0); +} + +// =========================================================================== +// Test 3: theta=0 vs theta>0 give different predictions on nonlinear data. +// +// Use a logistic-map library and query with a point BETWEEN library rows +// (not an exact member) so the distance weights actually matter. +// =========================================================================== +TEST(SMapProjection, ThetaZeroAndPositiveDifferOnNonlinearData) { + auto series = LogisticMap(0.4, 3.9, 50); + auto lib = MakeLib(series, 3); + + SMapProjection smap_global; + smap_global.library = &lib; + smap_global.embedding_dimension = 3; + smap_global.theta = 0.0; + + SMapProjection smap_local; + smap_local.library = &lib; + smap_local.embedding_dimension = 3; + smap_local.theta = 8.0; // strong localization + + // Midpoint between library rows 10 and 11 → not an exact library member + fims::Vector q(3); + for (size_t j = 0; j < 3; ++j) { + q[j] = 0.5 * (lib.at(10, j) + lib.at(11, j)); + } + + double pred_global = smap_global.predict_one(q); + double pred_local = smap_local.predict_one(q); + + // Global (theta=0) and strongly local (theta=8) should differ on chaotic data + EXPECT_NE(pred_global, pred_local); +} + +// =========================================================================== +// Test 4: Gaussian kernel vs exponential kernel give different predictions. +// Validates that the kernel field is wired through and changes the result. +// =========================================================================== +TEST(SMapProjection, GaussianAndExponentialKernelsDiffer) { + auto series = LogisticMap(0.4, 3.9, 50); + auto lib = MakeLib(series, 3); + + SMapProjection smap_exp; + smap_exp.library = &lib; + smap_exp.embedding_dimension = 3; + smap_exp.theta = 3.0; + smap_exp.kernel = SMapKernel::kExponential; + + SMapProjection smap_gauss; + smap_gauss.library = &lib; + smap_gauss.embedding_dimension = 3; + smap_gauss.theta = 3.0; + smap_gauss.kernel = SMapKernel::kGaussian; + + // Midpoint query — not a library row + fims::Vector q(3); + for (size_t j = 0; j < 3; ++j) { + q[j] = 0.5 * (lib.at(5, j) + lib.at(6, j)); + } + + double pred_exp = smap_exp.predict_one(q); + double pred_gauss = smap_gauss.predict_one(q); + + // The two kernels apply different decay functions → results should differ + EXPECT_NE(pred_exp, pred_gauss); +} + +// =========================================================================== +// Test 5: predict() fills the predictions vector with the correct size. +// =========================================================================== +TEST(SMapProjection, PredictFillsOutputVectorCorrectSize) { + auto series = LogisticMap(0.3, 3.7, 30); + auto lib = MakeLib(series, 2); + auto test = MakeLib(series, 2); + + SMapProjection smap; + smap.library = &lib; + smap.embedding_dimension = 2; + smap.theta = 1.0; + + smap.predict(test); + EXPECT_EQ(smap.predictions.size(), test.n_rows); +} + +// =========================================================================== +// Test 6: Throws std::runtime_error when library is null. +// =========================================================================== +TEST(SMapProjection, ThrowsOnNullLibrary) { + SMapProjection smap; + smap.embedding_dimension = 2; + smap.theta = 1.0; + + fims::Vector q = {1.0, 2.0}; + EXPECT_THROW(smap.predict_one(q), std::runtime_error); +} + +// =========================================================================== +// Test 7: Throws std::invalid_argument when query dimension mismatches. +// =========================================================================== +TEST(SMapProjection, ThrowsOnQueryDimensionMismatch) { + auto series = LogisticMap(0.4, 3.9, 20); + auto lib = MakeLib(series, 2); + + SMapProjection smap; + smap.library = &lib; + smap.embedding_dimension = 2; + smap.theta = 1.0; + + fims::Vector q = {1.0}; // dimension 1, needs 2 + EXPECT_THROW(smap.predict_one(q), std::invalid_argument); +} + +// =========================================================================== +// Test 8: Throws when library has fewer than E+1 rows. +// =========================================================================== +TEST(SMapProjection, ThrowsWhenLibraryTooSmall) { + // E=3, h=1 → need series.size() >= lag_span + h + 1 = 2 + 1 + 1 = 4 rows. + // Series of length 5 gives n_rows = 5 - 2 - 1 = 2 < E+1 = 4. + fims::Vector series = {1.0, 2.0, 3.0, 4.0, 5.0}; + auto lib = MakeLib(series, 3); + + SMapProjection smap; + smap.library = &lib; + smap.embedding_dimension = 3; + smap.theta = 1.0; + + fims::Vector q = {5.0, 4.0, 3.0}; + EXPECT_THROW(smap.predict_one(q), std::runtime_error); +} + +// =========================================================================== +// Test 9: predict() throws when test embedding dimension mismatches. +// =========================================================================== +TEST(SMapProjection, PredictThrowsOnTestDimensionMismatch) { + auto series = LogisticMap(0.4, 3.9, 20); + auto lib = MakeLib(series, 2); + auto test = MakeLib(series, 3); // wrong E + + SMapProjection smap; + smap.library = &lib; + smap.embedding_dimension = 2; + smap.theta = 1.0; + + EXPECT_THROW(smap.predict(test), std::invalid_argument); +} + +// =========================================================================== +// Test 10: Negative theta throws std::invalid_argument (from SMapWeights). +// =========================================================================== +TEST(SMapProjection, ThrowsOnNegativeTheta) { + auto series = LogisticMap(0.4, 3.9, 20); + auto lib = MakeLib(series, 2); + + SMapProjection smap; + smap.library = &lib; + smap.embedding_dimension = 2; + smap.theta = -1.0; + + fims::Vector q = {0.5, 0.4}; + EXPECT_THROW(smap.predict_one(q), std::invalid_argument); +} + +// =========================================================================== +// Test 11: Logistic map — predictions stay within a plausible range. +// Logistic map lives in (0,1); predictions should be close. +// =========================================================================== +TEST(SMapProjection, LogisticMapPredictionInRange) { + auto series = LogisticMap(0.2, 3.9, 35); + auto lib = MakeLib(series, 3); + + SMapProjection smap; + smap.library = &lib; + smap.embedding_dimension = 3; + smap.theta = 2.0; + + size_t last = lib.n_rows - 1; + fims::Vector q(3); + for (size_t j = 0; j < 3; ++j) q[j] = lib.at(last, j); + + double pred = smap.predict_one(q); + EXPECT_GT(pred, -0.5); + EXPECT_LT(pred, 1.5); +} + +// =========================================================================== +// Test 12: forecast_horizon=2 gives a different target than h=1. +// Validates that the forecast_horizon parameter propagates correctly. +// =========================================================================== +TEST(SMapProjection, ForecastHorizonTwoGivesDifferentTargetThanOne) { + auto series = LogisticMap(0.4, 3.9, 30); + auto lib_h1 = MakeLib(series, 2, 1, 1); // h=1 (default) + auto lib_h2 = MakeLib(series, 2, 1, 2); // h=2 + + // Sanity: targets for the same embedded row differ between h=1 and h=2 + // target_h1[row] = series[target_index + 1] + // target_h2[row] = series[target_index + 2] + EXPECT_NE(*lib_h1.target_values[0], *lib_h2.target_values[0]); +} + +// =========================================================================== +// Test 13: Extreme theta (ill-conditioned matrix) handled by Eigen::LDLT. +// Steve Munch / Andrea review: verify LDLT least-squares fallback. +// =========================================================================== +TEST(SMapProjection, IllConditionedExtremeThetaHandledByLDLT) { + auto series = LogisticMap(0.2, 3.9, 35); + auto lib = MakeLib(series, 3); + + SMapProjection smap; + smap.library = &lib; + smap.embedding_dimension = 3; + smap.theta = 1000.0; // Extremely large theta + + size_t last = lib.n_rows - 1; + fims::Vector q(3); + for (size_t j = 0; j < 3; ++j) q[j] = lib.at(last, j); + + EXPECT_NO_THROW({ + double pred = smap.predict_one(q); + EXPECT_FALSE(std::isnan(pred)); + EXPECT_FALSE(std::isinf(pred)); + }); +} + +} // namespace + diff --git a/tests/testthat/test-edm-fimsframe.R b/tests/testthat/test-edm-fimsframe.R index cd4123da1..be3c7f983 100644 --- a/tests/testthat/test-edm-fimsframe.R +++ b/tests/testthat/test-edm-fimsframe.R @@ -160,12 +160,12 @@ test_that("create_edm_embedding() stores uncertainty fields when uncertainty_nam # Grab the existing survey1 index rows and fabricate a matching sd series base_data <- get_data(data_4_model) |> dplyr::filter(.data[["type"]] == "index", - .data[["name"]] == "survey1") |> + .data[["fleet"]] == "survey1") |> dplyr::arrange(.data[["timing"]]) sd_rows <- base_data |> dplyr::mutate( - name = "survey1_sd", + fleet = "survey1_sd", value = 0.1 * .data[["value"]] ) diff --git a/tests/testthat/test-rcpp-edm.R b/tests/testthat/test-rcpp-edm.R index aea8bc43a..be52ff793 100644 --- a/tests/testthat/test-rcpp-edm.R +++ b/tests/testthat/test-rcpp-edm.R @@ -38,19 +38,20 @@ test_that("rcpp edm works with correct inputs", { #' @description Test that time_lag is set correctly after construction. expect_equal(de$time_lag, tau) #' @description Test that n_rows is set correctly. - expect_equal(de$n_rows, 3) + expect_equal(de$n_rows, 2) #' @description Test that n_cols is set correctly. expect_equal(de$n_cols, 3) #' @description Test target values mapping (actual x_t values, not indices). - expect_equal(de$target_values$toRVector(), c(30.0, 40.0, 50.0)) + expect_equal(de$target_values$get_values(), c(40.0, 50.0)) #' @description Test retrieving elements using at() method. expect_equal(de$at(0, 0), 30.0) expect_equal(de$at(0, 1), 20.0) expect_equal(de$at(0, 2), 10.0) expect_equal(de$at(1, 0), 40.0) - expect_equal(de$at(2, 2), 30.0) + expect_equal(de$at(1, 1), 30.0) + expect_equal(de$at(1, 2), 20.0) clear() }) @@ -87,10 +88,9 @@ test_that("rcpp edm construct_drop_missing works correctly", { de2$construct_drop_missing(series2, E, tau, missing_val) #' @description Test that construct_drop_missing keeps valid windows and maps correct target values. - expect_equal(de2$n_rows, 2) - expect_equal(de2$target_values$toRVector(), c(30.0, 70.0)) + expect_equal(de2$n_rows, 1) + expect_equal(de2$target_values$get_values(), -999.0) expect_equal(de2$at(0, 0), 30.0) - expect_equal(de2$at(1, 2), 50.0) clear() }) @@ -122,27 +122,24 @@ test_that("rcpp edm propagates uncertainty vectors when provided", { # Shape is identical to the value embedding #' @description Test that n_rows matches the value embedding when uncertainty is provided. - expect_equal(de$n_rows, 3) + expect_equal(de$n_rows, 2) #' @description Test that n_cols matches the value embedding when uncertainty is provided. expect_equal(de$n_cols, 3) # target_uncertainty mirrors target_values in the uncertainty space: - # Row 0 target index 2 -> sigma_2 = 0.3 - # Row 1 target index 3 -> sigma_3 = 0.4 - # Row 2 target index 4 -> sigma_4 = 0.5 + # Row 0 target index 2 -> sigma_3 = 0.4 + # Row 1 target index 3 -> sigma_4 = 0.5 #' @description Test that target_uncertainty holds the correct sigma_t values per row. - expect_equal(de$target_uncertainty$toRVector(), c(0.3, 0.4, 0.5)) + expect_equal(de$target_uncertainty$get_values(), c(0.4, 0.5)) # embedded_uncertainty row-major layout mirrors embedded_values: # Row 0: [sigma_2, sigma_1, sigma_0] = [0.3, 0.2, 0.1] # Row 1: [sigma_3, sigma_2, sigma_1] = [0.4, 0.3, 0.2] - # Row 2: [sigma_4, sigma_3, sigma_2] = [0.5, 0.4, 0.3] #' @description Test that embedded_uncertainty is laid out row-major matching embedded_values. expect_equal( - de$embedded_uncertainty$toRVector(), + de$embedded_uncertainty$get_values(), c(0.3, 0.2, 0.1, - 0.4, 0.3, 0.2, - 0.5, 0.4, 0.3) + 0.4, 0.3, 0.2) ) clear() @@ -155,9 +152,9 @@ test_that("rcpp edm uncertainty fields are empty when not provided", { de$construct(series, 3L, 1L) #' @description Test that embedded_uncertainty is empty when no uncertainty series is passed. - expect_equal(length(de$embedded_uncertainty$toRVector()), 0) + expect_equal(length(de$embedded_uncertainty$get_values()), 0) #' @description Test that target_uncertainty is empty when no uncertainty series is passed. - expect_equal(length(de$target_uncertainty$toRVector()), 0) + expect_equal(length(de$target_uncertainty$get_values()), 0) clear() }) @@ -175,12 +172,378 @@ test_that("rcpp edm construct_drop_missing propagates uncertainty correctly", { de$construct_drop_missing_with_uncertainty(series2, E, tau, missing_val, uncertainty) #' @description Test that construct_drop_missing retains only valid-row uncertainties. - expect_equal(de$n_rows, 2) + expect_equal(de$n_rows, 1) - # Row 0 target index 2 -> sigma_2 = 0.3 - # Row 1 target index 6 -> sigma_6 = 0.7 + # Row 0 target index 2 -> sigma_3 = 0.4 #' @description Test that target_uncertainty holds sigma_t for the retained rows. - expect_equal(de$target_uncertainty$toRVector(), c(0.3, 0.7)) + expect_equal(de$target_uncertainty$get_values(), c(0.4)) + + clear() +}) + +# rcpp edm SimplexProjection ---- +## Setup ---- +# Shared helper: build a DelayEmbeddingInterface from a numeric vector. +.make_de <- function(series, E = 3L, tau = 1L) { + de <- methods::new(DelayEmbedding) + de$construct(series, E, tau) + de +} + +## IO correctness ---- +test_that("SimplexProjection returns correct output length and finite values", { + series <- seq(0.1, 2.0, by = 0.1) # 20 evenly-spaced points + E <- 3L + lib_de <- .make_de(series[1:15], E) + test_de <- .make_de(series[13:20], E) + + sp <- methods::new(SimplexProjection) + sp$embedding_dimension <- E + sp$n_neighbors <- E + 1L + preds <- sp$predict(lib_de$get_id(), test_de$get_id()) + + #' @description Test that SimplexProjection::predict returns a numeric vector. + expect_true(is.numeric(preds)) + #' @description Test that SimplexProjection output length equals the number of test rows. + expect_equal(length(preds), test_de$n_rows) + #' @description Test that all SimplexProjection predictions are finite. + expect_true(all(is.finite(preds))) + + clear() +}) + +test_that("SimplexProjection matches pure-R squared-distance implementation", { + series <- seq(0.1, 2.0, by = 0.1) + E <- 2L; k <- E + 1L + lib_de <- .make_de(series[1:15], E) + test_de <- .make_de(series[13:20], E) + + sp <- methods::new(SimplexProjection) + sp$embedding_dimension <- E + sp$n_neighbors <- k + fims_pred <- sp$predict(lib_de$get_id(), test_de$get_id()) + + lib_mat <- matrix(lib_de$embedded_values$get_values(), ncol = E, byrow = TRUE) + lib_targets <- lib_de$target_values$get_values() + test_mat <- matrix(test_de$embedded_values$get_values(), ncol = E, byrow = TRUE) + r_pred <- vapply(seq_len(nrow(test_mat)), function(i) { + sq_d <- rowSums(sweep(lib_mat, 2, test_mat[i, ], "-")^2) + idx <- order(sq_d)[seq_len(k)] + d_min <- min(sq_d[idx]) + w <- exp(-sq_d[idx] / (d_min + 1e-12)) + w <- w / sum(w) + sum(w * lib_targets[idx]) + }, numeric(1)) + + #' @description Test that SimplexProjection matches a pure-R squared-distance implementation to within machine precision. + expect_equal(fims_pred, r_pred, tolerance = 1e-12) + + clear() +}) + +test_that("SimplexProjection predictions lie within the range of library targets", { + series <- c(0.1, 0.5, 0.9, 0.4, 0.7, 0.3, 0.8, 0.2, 0.6, 0.1, + 0.5, 0.9, 0.4, 0.7, 0.3, 0.8, 0.2, 0.6) + E <- 2L + lib_de <- .make_de(series[1:12], E) + test_de <- .make_de(series[10:18], E) + + sp <- methods::new(SimplexProjection) + sp$embedding_dimension <- E + sp$n_neighbors <- E + 1L + preds <- sp$predict(lib_de$get_id(), test_de$get_id()) + + lib_targets <- lib_de$target_values$get_values() + #' @description Test that SimplexProjection predictions stay within the library target range (weighted-average property). + expect_true(all(preds >= min(lib_targets) - 1e-10)) + expect_true(all(preds <= max(lib_targets) + 1e-10)) + + clear() +}) + +## Error handling ---- +test_that("SimplexProjection throws on invalid library id", { + series <- seq(0.1, 0.6, by = 0.1) + test_de <- .make_de(series, 2L) + + sp <- methods::new(SimplexProjection) + sp$embedding_dimension <- 2L + sp$n_neighbors <- 3L + + #' @description Test that SimplexProjection::predict throws when given a non-existent library id. + expect_error(sp$predict(9999L, test_de$get_id())) + + clear() +}) + +test_that("SimplexProjection throws on invalid test id", { + series <- seq(0.1, 0.6, by = 0.1) + lib_de <- .make_de(series, 2L) + + sp <- methods::new(SimplexProjection) + sp$embedding_dimension <- 2L + sp$n_neighbors <- 3L + + #' @description Test that SimplexProjection::predict throws when given a non-existent test id. + expect_error(sp$predict(lib_de$get_id(), 9999L)) + + clear() +}) + +# rcpp edm SMapProjection ---- +## IO correctness ---- +test_that("SMapProjection returns correct output length and finite values", { + # Chaotic-like series avoids collinear embedding columns in WLS design matrix + series <- c(0.4, 0.9, 0.3, 0.8, 0.5, 0.7, 0.2, 0.6, 0.4, 0.9, + 0.3, 0.8, 0.5, 0.7, 0.2, 0.6, 0.4, 0.9, 0.3, 0.8) + E <- 2L + lib_de <- .make_de(series[1:15], E) + test_de <- .make_de(series[13:20], E) + + sm <- methods::new(SMapProjection) + sm$embedding_dimension <- E + sm$theta <- 1.0 + sm$kernel <- "exponential" + preds <- sm$predict(lib_de$get_id(), test_de$get_id()) + + #' @description Test that SMapProjection::predict returns a numeric vector. + expect_true(is.numeric(preds)) + #' @description Test that SMapProjection output length equals the number of test rows. + expect_equal(length(preds), test_de$n_rows) + #' @description Test that all SMapProjection predictions are finite. + expect_true(all(is.finite(preds))) + + clear() +}) + +test_that("SMapProjection matches pure-R WLS squared-distance implementation", { + # Chaotic-like series avoids collinear embedding columns in WLS design matrix + series <- c(0.4, 0.9, 0.3, 0.8, 0.5, 0.7, 0.2, 0.6, 0.4, 0.9, + 0.3, 0.8, 0.5, 0.7, 0.2, 0.6, 0.4, 0.9, 0.3, 0.8) + E <- 2L; theta <- 1.5 + lib_de <- .make_de(series[1:15], E) + test_de <- .make_de(series[13:20], E) + + sm <- methods::new(SMapProjection) + sm$embedding_dimension <- E + sm$theta <- theta + sm$kernel <- "exponential" + fims_pred <- sm$predict(lib_de$get_id(), test_de$get_id()) + + lib_mat <- matrix(lib_de$embedded_values$get_values(), ncol = E, byrow = TRUE) + lib_targets <- lib_de$target_values$get_values() + test_mat <- matrix(test_de$embedded_values$get_values(), ncol = E, byrow = TRUE) + X <- cbind(1, lib_mat) + r_pred <- vapply(seq_len(nrow(test_mat)), function(i) { + sq_d <- rowSums(sweep(lib_mat, 2, test_mat[i, ], "-")^2) + d_mean <- mean(sq_d) + w <- exp(-theta * sq_d / (d_mean + 1e-12)) + coef <- solve(t(X) %*% diag(w) %*% X, t(X) %*% diag(w) %*% lib_targets) + sum(coef * c(1, test_mat[i, ])) + }, numeric(1)) + + #' @description Test that SMapProjection matches a pure-R WLS squared-distance implementation to within 1e-10. + expect_equal(fims_pred, r_pred, tolerance = 1e-10) + + clear() +}) + +test_that("SMapProjection gaussian kernel gives different predictions than exponential", { + # Chaotic-like series avoids collinear embedding columns in WLS design matrix + series <- c(0.4, 0.9, 0.3, 0.8, 0.5, 0.7, 0.2, 0.6, 0.4, 0.9, + 0.3, 0.8, 0.5, 0.7, 0.2, 0.6, 0.4, 0.9, 0.3, 0.8) + E <- 2L + lib_de <- .make_de(series[1:14], E) + test_de <- .make_de(series[12:20], E) + + sm_exp <- methods::new(SMapProjection) + sm_exp$embedding_dimension <- E + sm_exp$theta <- 1.0 + sm_exp$kernel <- "exponential" + pred_exp <- sm_exp$predict(lib_de$get_id(), test_de$get_id()) + + sm_gauss <- methods::new(SMapProjection) + sm_gauss$embedding_dimension <- E + sm_gauss$theta <- 1.0 + sm_gauss$kernel <- "gaussian" + pred_gauss <- sm_gauss$predict(lib_de$get_id(), test_de$get_id()) + + #' @description Test that gaussian and exponential SMap kernels give different predictions. + expect_false(isTRUE(all.equal(pred_exp, pred_gauss))) + + clear() +}) + +## Error handling ---- +test_that("SMapProjection throws on invalid library id", { + series <- seq(0.1, 1.0, by = 0.1) + test_de <- .make_de(series, 2L) + + sm <- methods::new(SMapProjection) + sm$embedding_dimension <- 2L + sm$theta <- 1.0 + + #' @description Test that SMapProjection::predict throws when given a non-existent library id. + expect_error(sm$predict(9999L, test_de$get_id())) + + clear() +}) + +# rcpp edm GPEdmProjection ---- +## IO correctness ---- +test_that("GPEdmProjection predict returns correct output length and finite values", { + series <- seq(0.1, 1.5, by = 0.1) # 15 points + E <- 2L + lib_de <- .make_de(series[1:10], E) + test_de <- .make_de(series[8:15], E) + + gp <- methods::new(GPEdmProjection) + gp$embedding_dimension <- E + gp$phi <- rep(0.5, E) + gp$sigma2 <- 1.0 + gp$ve <- 0.1 + preds <- gp$predict(lib_de$get_id(), test_de$get_id()) + + #' @description Test that GPEdmProjection::predict returns a numeric vector. + expect_true(is.numeric(preds)) + #' @description Test that GPEdmProjection output length equals the number of test rows. + expect_equal(length(preds), test_de$n_rows) + #' @description Test that all GPEdmProjection predictions are finite. + expect_true(all(is.finite(preds))) + + clear() +}) + +test_that("GPEdmProjection fit returns named list with phi, sigma2, ve", { + series <- c(0.4, 0.9, 0.3, 0.8, 0.5, 0.7, 0.2, 0.6, 0.4, 0.9, + 0.3, 0.8, 0.5, 0.7, 0.2) + E <- 2L + lib_de <- .make_de(series, E) + + gp <- methods::new(GPEdmProjection) + gp$embedding_dimension <- E + gp$phi <- rep(0.1, E) + gp$sigma2 <- 1.0 + gp$ve <- 0.1 + result <- gp$fit(lib_de$get_id()) + + #' @description Test that GPEdmProjection::fit returns a list. + expect_true(is.list(result)) + #' @description Test that the fit result contains the phi hyperparameter vector. + expect_true("phi" %in% names(result)) + #' @description Test that the fit result contains the sigma2 hyperparameter. + expect_true("sigma2" %in% names(result)) + #' @description Test that the fit result contains the ve hyperparameter. + expect_true("ve" %in% names(result)) + #' @description Test that fitted phi has length equal to embedding_dimension. + expect_equal(length(result$phi), E) + #' @description Test that all fitted hyperparameters are finite. + expect_true(all(is.finite(result$phi))) + expect_true(is.finite(result$sigma2)) + expect_true(is.finite(result$ve)) + + clear() +}) + +test_that("GPEdmProjection fit updates interface fields with optimized hyperparameters", { + series <- c(0.4, 0.9, 0.3, 0.8, 0.5, 0.7, 0.2, 0.6, 0.4, 0.9, + 0.3, 0.8, 0.5, 0.7, 0.2) + E <- 2L + lib_de <- .make_de(series, E) + + gp <- methods::new(GPEdmProjection) + gp$embedding_dimension <- E + gp$phi <- rep(0.1, E) + gp$sigma2 <- 1.0 + gp$ve <- 0.1 + gp$fit(lib_de$get_id()) + + #' @description Test that GPEdmProjection::fit writes the optimized phi back to the interface object. + expect_equal(length(gp$phi), E) + expect_true(all(is.finite(gp$phi))) + #' @description Test that GPEdmProjection::fit writes the optimized sigma2 back to the interface object. + expect_true(is.finite(gp$sigma2)) + #' @description Test that GPEdmProjection::fit writes the optimized ve back to the interface object. + expect_true(is.finite(gp$ve)) + + clear() +}) + +test_that("GPEdmProjection fit followed by predict gives finite predictions", { + series <- c(0.4, 0.9, 0.3, 0.8, 0.5, 0.7, 0.2, 0.6, 0.4, 0.9, + 0.3, 0.8, 0.5, 0.7, 0.2, 0.6, 0.4, 0.9, 0.3, 0.8) + E <- 2L + lib_de <- .make_de(series[1:14], E) + test_de <- .make_de(series[12:20], E) + + gp <- methods::new(GPEdmProjection) + gp$embedding_dimension <- E + gp$phi <- rep(0.1, E) + gp$sigma2 <- 1.0 + gp$ve <- 0.1 + gp$fit(lib_de$get_id()) + preds <- gp$predict(lib_de$get_id(), test_de$get_id()) + + #' @description Test that fit-then-predict gives finite GP-EDM predictions. + expect_true(all(is.finite(preds))) + #' @description Test that fit-then-predict output length equals the number of test embedding rows. + expect_equal(length(preds), test_de$n_rows) + + clear() +}) + +test_that("GPEdmProjection predictions differ for small vs. large phi (kernel localization)", { + series <- seq(0.05, 1.0, by = 0.05) # 20 smooth points + E <- 2L + lib_de <- .make_de(series[1:14], E) + test_de <- .make_de(series[12:20], E) + + gp_small <- methods::new(GPEdmProjection) + gp_small$embedding_dimension <- E + gp_small$phi <- rep(0.01, E) + gp_small$sigma2 <- 1.0 + gp_small$ve <- 0.01 + pred_small <- gp_small$predict(lib_de$get_id(), test_de$get_id()) + + gp_large <- methods::new(GPEdmProjection) + gp_large$embedding_dimension <- E + gp_large$phi <- rep(100.0, E) + gp_large$sigma2 <- 1.0 + gp_large$ve <- 0.01 + pred_large <- gp_large$predict(lib_de$get_id(), test_de$get_id()) + + #' @description Test that GP-EDM predictions differ between small and large phi (kernel localization). + expect_false(isTRUE(all.equal(pred_small, pred_large, tolerance = 1e-4))) + + clear() +}) + +## Error handling ---- +test_that("GPEdmProjection predict throws on invalid library id", { + series <- seq(0.1, 1.0, by = 0.1) + test_de <- .make_de(series, 2L) + + gp <- methods::new(GPEdmProjection) + gp$embedding_dimension <- 2L + gp$phi <- rep(0.5, 2L) + gp$sigma2 <- 1.0 + gp$ve <- 0.1 + + #' @description Test that GPEdmProjection::predict throws when given a non-existent library id. + expect_error(gp$predict(9999L, test_de$get_id())) + + clear() +}) + +test_that("GPEdmProjection fit throws on invalid library id", { + gp <- methods::new(GPEdmProjection) + gp$embedding_dimension <- 2L + gp$phi <- rep(0.1, 2L) + gp$sigma2 <- 1.0 + gp$ve <- 0.1 + + #' @description Test that GPEdmProjection::fit throws when given a non-existent library id. + expect_error(gp$fit(9999L)) clear() })