draft rcpp wrapper for dmi

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)
 }

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<double> dmi(
         const std::vector<double>& vec,
-        const std::vector<size_t>& pred,
         const std::vector<size_t>& tau,
+        const std::vector<size_t>& pred,
         size_t k = 3,
         size_t alg = 0,
         double base = 2.0,

src/DMI.cpp

@@ -0,0 +1,83 @@
+#include <vector>
+#include <cmath>
+#include <limits>
+#include <string>
+#include <utility>
+#include <numeric>
+#include <algorithm>
+#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<double> tg = Rcpp::as<std::vector<double>>(target);
+    const size_t n_obs = tg.size();
+
+    // Convert prediction indices (R 1-based → C++ 0-based)
+    std::vector<size_t> pred_std = Rcpp::as<std::vector<size_t>>(pred);
+    for (auto& idx : pred_std) 
+    {
+        if (idx < 1 || idx > n_obs) 
+        {
+            Rcpp::stop("pred index %d out of bounds [1, %d]",
+                       static_cast<int>(idx),
+                       static_cast<int>(n_obs));
+        }
+        idx -= 1;
+    }
+
+    // Construct time delay step tau
+    std::vector<size_t> tau_std = Rcpp::as<std::vector<size_t>>(tau);
+    if (tau_std.empty()) {
+        Rcpp::stop("tau vector cannot be empty.");
+    }
+    size_t max_tau = static_cast<size_t>(*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<double> res = pc::dmi::dmi(
+            tg, tau_std, pred_std, 
+            static_cast<size_t>(std::abs(k)), 
+            static_cast<size_t>(std::abs(alg)), 
+            base, normalize,
+            static_cast<size_t>(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;
+}

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<size_t> 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<double> 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<size_t>(std::abs(k)), 

src/RcppExports.cpp

@@ -12,6 +12,23 @@ Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
 Rcpp::Rostream<false>& 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<Rcpp::List> nb, Rcpp::Nullable<int> 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},