diff --git a/R/RcppExports.R b/R/RcppExports.R index 7b3f4f2..7eab922 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,6 +1,10 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 +RcppDMI <- function(target, tau, pred, k = 3L, alg = 0L, base = 2.0, normalize = FALSE, threads = 1L) { + .Call(`_pc_RcppDMI`, target, tau, pred, k, alg, base, normalize, threads) +} + RcppFNN <- function(target, rt, eps, lib, pred, E, tau = 1L, style = 0L, dist_metric = "euclidean", k = 3L, threads = 1L, parallel_level = 0L, nb = NULL, nrows = NULL) { .Call(`_pc_RcppFNN`, target, rt, eps, lib, pred, E, tau, style, dist_metric, k, threads, parallel_level, nb, nrows) } diff --git a/inst/include/pc/dmi.hpp b/inst/include/pc/dmi.hpp index 1b614f6..f5b65a9 100644 --- a/inst/include/pc/dmi.hpp +++ b/inst/include/pc/dmi.hpp @@ -48,8 +48,8 @@ namespace dmi * * Parameters: * vec - Input numeric vector representing the ordered series - * pred - Indices defining the sample positions (from past to present) * tau - Vector of lag steps (non-negative integers) + * pred - Indices defining the sample positions (from past to present) * k - Number of nearest neighbors for KSG estimator (default: 3) * alg - Algorithm variant for KSG estimator (default: 0) * base - Logarithm base for mutual information (default: 2.0) @@ -63,8 +63,8 @@ namespace dmi */ inline std::vector dmi( const std::vector& vec, - const std::vector& pred, const std::vector& tau, + const std::vector& pred, size_t k = 3, size_t alg = 0, double base = 2.0, diff --git a/src/DMI.cpp b/src/DMI.cpp new file mode 100644 index 0000000..a7337c7 --- /dev/null +++ b/src/DMI.cpp @@ -0,0 +1,83 @@ +#include +#include +#include +#include +#include +#include +#include +#include "pc.h" + +// Wrapper function to perform delayed mutual information analysis +// [[Rcpp::export(rng = false)]] +Rcpp::NumericVector RcppDMI( + const Rcpp::NumericVector& target, + const Rcpp::NumericVector& tau, + const Rcpp::IntegerVector& pred, + int k = 3, + int alg = 0, + double base = 2.0, + bool normalize = false, + int threads = 1) +{ + // --- Input Conversion and Validation --- + std::vector tg = Rcpp::as>(target); + const size_t n_obs = tg.size(); + + // Convert prediction indices (R 1-based → C++ 0-based) + std::vector pred_std = Rcpp::as>(pred); + for (auto& idx : pred_std) + { + if (idx < 1 || idx > n_obs) + { + Rcpp::stop("pred index %d out of bounds [1, %d]", + static_cast(idx), + static_cast(n_obs)); + } + idx -= 1; + } + + // Construct time delay step tau + std::vector tau_std = Rcpp::as>(tau); + if (tau_std.empty()) { + Rcpp::stop("tau vector cannot be empty."); + } + size_t max_tau = static_cast(*std::max_element(tau_std.begin(), tau_std.end())); + + // ---- sort predict indices ---- + pred_std.erase( + std::remove_if(pred_std.begin(), pred_std.end(), + [&](size_t idx){ return idx < max_tau; }), + pred_std.end() + ); + + std::sort(pred_std.begin(), pred_std.end()); + pred_std.erase( + std::unique(pred_std.begin(), pred_std.end()), + pred_std.end() + ); + + // ---- filter pred (remove NaN in target) ---- + pred_std.erase( + std::remove_if(pred_std.begin(), pred_std.end(), + [&](size_t idx){ return std::isnan(tg[idx]); }), + pred_std.end() + ); + + // --- Perform Delay Mutual Information Analysis --- + std::vector res = pc::dmi::dmi( + tg, tau_std, pred_std, + static_cast(std::abs(k)), + static_cast(std::abs(alg)), + base, normalize, + static_cast(std::abs(threads))); + + // Convert the result back to Rcpp::NumericVector and set names as "tau:1", "tau:2", ..., "tau:n" + Rcpp::NumericVector result = Rcpp::wrap(res); + Rcpp::CharacterVector resnames(result.size()); + for (int i = 0; i < result.size(); ++i) { + resnames[i] = "tau:" + std::to_string(tau_std[i]); + } + result.names() = resnames; + + return result; +} diff --git a/src/FNN.cpp b/src/FNN.cpp index 13c2738..230ffc2 100644 --- a/src/FNN.cpp +++ b/src/FNN.cpp @@ -160,28 +160,18 @@ Rcpp::NumericVector RcppFNN( pred_std.end() ); - // ---- filter lib/pred (remove NaN in target/source) ---- - size_t write = 0; - for (size_t i = 0; i < lib_std.size(); ++i) - { - size_t idx = lib_std[i]; - if (!std::isnan(tg[idx])) - { - lib_std[write++] = idx; - } - } - lib_std.resize(write); + // ---- filter lib/pred (remove NaN in target) ---- + lib_std.erase( + std::remove_if(lib_std.begin(), lib_std.end(), + [&](size_t idx){ return std::isnan(tg[idx]); }), + lib_std.end() + ); - write = 0; - for (size_t i = 0; i < pred_std.size(); ++i) - { - size_t idx = pred_std[i]; - if (!std::isnan(tg[idx])) - { - pred_std[write++] = idx; - } - } - pred_std.resize(write); + pred_std.erase( + std::remove_if(pred_std.begin(), pred_std.end(), + [&](size_t idx){ return std::isnan(tg[idx]); }), + pred_std.end() + ); // --- Prepare for data slicing --- std::vector selected_indices; @@ -199,7 +189,7 @@ Rcpp::NumericVector RcppFNN( // --- Check if full set is used --- bool use_subset = (selected_indices.size() < Mx.size()); - // --- Perform Pattern Causality Analysis --- + // --- Perform FNN Analysis --- std::vector res; if (!use_subset) @@ -244,7 +234,7 @@ Rcpp::NumericVector RcppFNN( pred_std[i] = index_map[pred_std[i]]; } - // --- Run patcaus on subset --- + // --- Run fnn on subset --- res = pc::fnn::fnn( Mx_sub, lib_std, pred_std, rt_std, eps_std, dist_metric, static_cast(std::abs(k)), diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 902bf09..8fcfe42 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -12,6 +12,23 @@ Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif +// RcppDMI +Rcpp::NumericVector RcppDMI(const Rcpp::NumericVector& target, const Rcpp::NumericVector& tau, const Rcpp::IntegerVector& pred, int k, int alg, double base, bool normalize, int threads); +RcppExport SEXP _pc_RcppDMI(SEXP targetSEXP, SEXP tauSEXP, SEXP predSEXP, SEXP kSEXP, SEXP algSEXP, SEXP baseSEXP, SEXP normalizeSEXP, SEXP threadsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type target(targetSEXP); + Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type tau(tauSEXP); + Rcpp::traits::input_parameter< const Rcpp::IntegerVector& >::type pred(predSEXP); + Rcpp::traits::input_parameter< int >::type k(kSEXP); + Rcpp::traits::input_parameter< int >::type alg(algSEXP); + Rcpp::traits::input_parameter< double >::type base(baseSEXP); + Rcpp::traits::input_parameter< bool >::type normalize(normalizeSEXP); + Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); + rcpp_result_gen = Rcpp::wrap(RcppDMI(target, tau, pred, k, alg, base, normalize, threads)); + return rcpp_result_gen; +END_RCPP +} // RcppFNN Rcpp::NumericVector RcppFNN(const Rcpp::NumericVector& target, const Rcpp::NumericVector& rt, const Rcpp::NumericVector& eps, const Rcpp::IntegerVector& lib, const Rcpp::IntegerVector& pred, const Rcpp::IntegerVector& E, int tau, int style, const std::string& dist_metric, int k, int threads, int parallel_level, Rcpp::Nullable nb, Rcpp::Nullable nrows); RcppExport SEXP _pc_RcppFNN(SEXP targetSEXP, SEXP rtSEXP, SEXP epsSEXP, SEXP libSEXP, SEXP predSEXP, SEXP ESEXP, SEXP tauSEXP, SEXP styleSEXP, SEXP dist_metricSEXP, SEXP kSEXP, SEXP threadsSEXP, SEXP parallel_levelSEXP, SEXP nbSEXP, SEXP nrowsSEXP) { @@ -120,6 +137,7 @@ END_RCPP } static const R_CallMethodDef CallEntries[] = { + {"_pc_RcppDMI", (DL_FUNC) &_pc_RcppDMI, 8}, {"_pc_RcppFNN", (DL_FUNC) &_pc_RcppFNN, 14}, {"_pc_RcppPC", (DL_FUNC) &_pc_RcppPC, 16}, {"_pc_RcppPCboot", (DL_FUNC) &_pc_RcppPCboot, 22},