Merge pull request #17 from stscl/dev

update `gccm`
stscl · Dec 10, 2024 · 2310f0a · 2310f0a
2 parents f51612e + a470842
commit 2310f0a
Show file tree

Hide file tree

Showing 15 changed files with 3,403 additions and 440 deletions.
diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -17,22 +17,6 @@ RcppGenEmbeddings <- function(vec, nbmat, E) {
     .Call(`_spEDM_RcppGenEmbeddings`, vec, nbmat, E)
 }
 
-RcppDist <- function(matrix) {
-    .Call(`_spEDM_RcppDist`, matrix)
-}
-
-RcppClosestIndices <- function(distmat, libsize) {
-    .Call(`_spEDM_RcppClosestIndices`, distmat, libsize)
-}
-
-RcppCCMWeight <- function(distmat, closestIndices, libsize) {
-    .Call(`_spEDM_RcppCCMWeight`, distmat, closestIndices, libsize)
-}
-
-RcppSimplexProjection <- function(y, x, nbmat, libsize, E) {
-    .Call(`_spEDM_RcppSimplexProjection`, y, x, nbmat, libsize, E)
-}
-
 RcppGCCMLattice <- function(y, x, nbmat, libsizes, E) {
     .Call(`_spEDM_RcppGCCMLattice`, y, x, nbmat, libsizes, E)
 }

diff --git a/R/gccm.R b/R/gccm.R
@@ -6,16 +6,16 @@
 #' @param coords (optional) The coordinates matrix.
 #' @param data (optional) The observation data.
 #' @param libsizes (optional) A vector of library sizes to use.
-#' @param E (optional) The number of dimensions for the attractor reconstruction.
+#' @param E (optional) The dimensions of the embedding.
 #' @param trend_rm (optional) Variables that need to remove linear trends.
 #' @param ... (optional) Other parameters passed to `sdsfun::spdep_nb()`.
 #'
 #' @return A `tibble`.
 #' @export
 #'
 #' @examples
-#' columbus <- sf::read_sf(system.file("shapes/columbus.gpkg", package="spData")[1],
-#'                         quiet=TRUE)
+#' columbus = sf::read_sf(system.file("shapes/columbus.gpkg", package="spData")[1],
+#'                        quiet=TRUE)
 #' gccm("CRIME","HOVAL", data = columbus)
 gccm = \(cause, effect, nbmat = NULL, coords = NULL, data = NULL, libsizes = NULL,
          E = 3, trend_rm = c("effect", "cause", "none", "all"), ...) {
@@ -45,7 +45,7 @@ gccm = \(cause, effect, nbmat = NULL, coords = NULL, data = NULL, libsizes = NUL
   }
 
   if (length(cause) != nrow(nbmat)) stop("Incompatible Data Dimensions!")
-  if (is.null(libsizes)) libsizes = floor(seq(1,length(cause),length.out = 15))
+  if (is.null(libsizes)) libsizes = floor(seq(E + 2,length(cause),length.out = 15))
 
   trend_rm = match.arg(trend_rm)
   switch(trend_rm,
@@ -62,12 +62,12 @@ gccm = \(cause, effect, nbmat = NULL, coords = NULL, data = NULL, libsizes = NUL
 
 
   x_xmap_y = RcppGCCMLattice(cause,effect,nbmat,libsizes,E)
-  colnames(x_xmap_y) = c("lib_sizes","x_xmap_y","x_xmap_y_sig",
+  colnames(x_xmap_y) = c("lib_sizes","x_xmap_y_mean","x_xmap_y_sig",
                          "x_xmap_y_upper","x_xmap_y_lower")
   x_xmap_y = tibble::as_tibble(x_xmap_y)
 
   y_xmap_x = RcppGCCMLattice(effect,cause,nbmat,libsizes,E)
-  colnames(y_xmap_x) = c("lib_sizes","y_xmap_x","y_xmap_x_sig",
+  colnames(y_xmap_x) = c("lib_sizes","y_xmap_x_mean","y_xmap_x_sig",
                          "y_xmap_x_upper","y_xmap_x_lower")
   y_xmap_x = tibble::as_tibble(y_xmap_x)
 

diff --git a/inst/extdata/popdensity.csv b/inst/extdata/popdensity.csv
diff --git a/man/figures/gccm/gccm_plot-1.png b/man/figures/gccm/gccm_plot-1.png
diff --git a/man/gccm.Rd b/man/gccm.Rd
diff --git a/src/CppGCCM.cpp b/src/CppGCCM.cpp
diff --git a/src/CppGCCM.h b/src/CppGCCM.h
diff --git a/src/CppUtils.cpp b/src/CppUtils.cpp
@@ -109,131 +109,3 @@ std::vector<std::vector<double>> GenEmbeddings(const std::vector<double>& vec,
 
   return embeddings;
 }
-
-// Function to calculate the pairwise absolute difference mean matrix
-std::vector<std::vector<double>> CppDist(const std::vector<std::vector<double>>& matrix) {
-  size_t n = matrix.size();
-  std::vector<std::vector<double>> result(n, std::vector<double>(n, std::numeric_limits<double>::quiet_NaN()));
-
-  for (size_t i = 0; i < n; ++i) {
-    for (size_t j = i; j < n; ++j) {
-      if (i == j) {
-        result[i][j] = 0.0; // Diagonal elements are zero
-      } else {
-        std::vector<double> abs_diff = CppAbs(matrix[i], matrix[j]);
-        result[i][j] = CppMean(abs_diff, true);
-        result[j][i] = result[i][j]; // Symmetric matrix
-      }
-    }
-  }
-
-  return result;
-}
-
-// Function to find the indices of the libsize + 1 closest elements for each row in distmat
-std::vector<std::vector<int>> CppClosestIndices(const std::vector<std::vector<double>>& distmat,
-                                                int libsize) {
-  size_t n = distmat.size();
-  std::vector<std::vector<int>> closestIndices(n, std::vector<int>(libsize + 1));
-
-  for (size_t i = 0; i < n; ++i) {
-    // Create a vector to store the indices and distances
-    std::vector<std::pair<double, int>> distances;
-
-    for (size_t j = 0; j < n; ++j) {
-      if (i != j) { // Exclude the diagonal element
-        distances.push_back(std::make_pair(distmat[i][j], j));
-      }
-    }
-
-    // Sort the distances vector by the distance (first element of the pair)
-    std::sort(distances.begin(), distances.end(),
-              [](const std::pair<double, int>& a,
-                 const std::pair<double, int>& b) {
-      return a.first < b.first;
-    });
-
-    // Select the libsize + 1 closest indices
-    for (int k = 0; k < libsize + 1; ++k) {
-      closestIndices[i][k] = distances[k].second;
-    }
-  }
-
-  return closestIndices;
-}
-
-// Function to calculate the CCM weights
-std::vector<std::vector<double>> CCMWeight(const std::vector<std::vector<double>>& distmat,
-                                           const std::vector<std::vector<int>>& closestIndices,
-                                           int libsize) {
-  size_t n = distmat.size();
-  std::vector<std::vector<double>> result(n, std::vector<double>(libsize + 1));
-
-  for (size_t i = 0; i < n; ++i) {
-    std::vector<double> sdist;
-
-    // Extract the corresponding distances from distmat
-    for (int j = 0; j < libsize + 1; ++j) {
-      int index = closestIndices[i][j];
-      sdist.push_back(distmat[i][index]);
-    }
-
-    // Normalize sdist by dividing each element by the first element,
-    // taking the negative and applying the exponential function
-    double firstElement = sdist[0];
-    for (size_t j = 0; j < sdist.size(); ++j) {
-      sdist[j] = std::exp(-sdist[j] / firstElement);
-    }
-
-    // Normalize the resulting vector using CppSumNormalize with NA_rm set to true
-    std::vector<double> normalizedWeights = CppSumNormalize(sdist, true);
-
-    // Assign the normalized weights to the result
-    for (int j = 0; j < libsize + 1; ++j) {
-      result[i][j] = normalizedWeights[j];
-    }
-  }
-
-  return result;
-}
-
-// Function to perform Simplex Projection
-std::vector<double> SimplexProjection(const std::vector<double>& y,
-                                      const std::vector<double>& x,
-                                      const std::vector<std::vector<int>>& nbmat,
-                                      int libsize,
-                                      int E) {
-  // Generate embeddings for x
-  std::vector<std::vector<double>> singleE = GenEmbeddings(x, nbmat, E);
-
-  // Calculate the distance matrix for singleE
-  std::vector<std::vector<double>> singleDist = CppDist(singleE);
-
-  // Find the closest indices for each row in singleDist
-  std::vector<std::vector<int>> singleCI = CppClosestIndices(singleDist, libsize);
-
-  // Calculate the weights for each row in singleCI
-  std::vector<std::vector<double>> singleW = CCMWeight(singleDist, singleCI, libsize);
-
-  // Construct singleY by extracting values from y based on singleCI
-  size_t n = y.size();
-  std::vector<std::vector<double>> singleY(n, std::vector<double>(libsize + 1));
-  for (size_t i = 0; i < n; ++i) {
-    for (int j = 0; j < libsize + 1; ++j) {
-      int index = singleCI[i][j];
-      singleY[i][j] = y[index];
-    }
-  }
-
-  // Calculate the result by multiplying singleY and singleW and summing the results
-  std::vector<double> result(n);
-  for (size_t i = 0; i < n; ++i) {
-    std::vector<double> sumvec;
-    for (int j = 0; j < libsize + 1; ++j) {
-      sumvec.push_back(singleY[i][j] * singleW[i][j]);
-    }
-    result[i] = CppSum(sumvec, true);
-  }
-
-  return result;
-}
diff --git a/src/CppUtils.h b/src/CppUtils.h
@@ -16,19 +16,4 @@ std::vector<std::vector<double>> GenEmbeddings(const std::vector<double>& vec,
                                                const std::vector<std::vector<int>>& nbmat,
                                                int E);
 
-std::vector<std::vector<double>> CppDist(const std::vector<std::vector<double>>& matrix);
-
-std::vector<std::vector<int>> CppClosestIndices(const std::vector<std::vector<double>>& distmat,
-                                                int libsize);
-
-std::vector<std::vector<double>> CCMWeight(const std::vector<std::vector<double>>& distmat,
-                                           const std::vector<std::vector<int>>& closestIndices,
-                                           int libsize);
-
-std::vector<double> SimplexProjection(const std::vector<double>& y,
-                                      const std::vector<double>& x,
-                                      const std::vector<std::vector<int>>& nbmat,
-                                      int libsize,
-                                      int E);
-
 #endif // CppUtils_H