diff --git a/DESCRIPTION b/DESCRIPTION index efef096..f35d104 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -21,8 +21,10 @@ LinkingTo: Imports: dplyr, sdsfun (>= 0.6.0), + terra, tibble Suggests: + ggplot2, knitr, Rcpp, RcppThread, diff --git a/R/RcppExports.R b/R/RcppExports.R index 1fcb1ec..a5d0c5c 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -9,6 +9,10 @@ RcppGenGridEmbeddings <- function(mat, E) { .Call(`_spEDM_RcppGenGridEmbeddings`, mat, E) } +RcppGCCM4Grid <- function(xMatrix, yMatrix, lib_sizes, pred, E) { + .Call(`_spEDM_RcppGCCM4Grid`, xMatrix, yMatrix, lib_sizes, pred, E) +} + RcppConfidence <- function(r, n, level = 0.05) { .Call(`_spEDM_RcppConfidence`, r, n, level) } diff --git a/src/CppGridUtils.cpp b/src/CppGridUtils.cpp index 8ed2ff9..a2f454d 100644 --- a/src/CppGridUtils.cpp +++ b/src/CppGridUtils.cpp @@ -3,105 +3,49 @@ #include #include -/** - * Expands a matrix by adding rows and columns of NaN values based on the lag number. - * - * @param mat The input matrix represented as a 2D vector of doubles. - * @param lagNum The number of times the matrix should be expanded. Must be a non-negative integer. - * @return The expanded matrix. - */ -std::vector> ExpandMat(std::vector> mat, - int lagNum) { - // If lagNum is negative, return the original matrix as no expansion is needed. - if (lagNum < 0) { - return mat; - } - - // If lagNum is greater than 1, recursively call ExpandMat to expand the matrix. - if (lagNum > 1) { - mat = ExpandMat(mat, lagNum - 1); - } - - // Get the number of columns and rows in the matrix. - int numCols = mat[0].size(); - // int numRows = mat.size(); - - // Add a row of NaN values at the top and bottom of the matrix. - std::vector nanRow(numCols, std::numeric_limits::quiet_NaN()); - mat.insert(mat.begin(), nanRow); // Add at the top - mat.push_back(nanRow); // Add at the bottom - - // Add a column of NaN values at the left and right of the matrix. - for (auto& row : mat) { - row.insert(row.begin(), std::numeric_limits::quiet_NaN()); // Add at the left - row.push_back(std::numeric_limits::quiet_NaN()); // Add at the right - } - - return mat; -} - -/** - * Generates a 2D matrix of lagged variables based on the input matrix and lag number. - * - * @param mat The input matrix represented as a 2D vector of doubles. - * @param lagNum The lag number used to generate the lagged variables. - * @return A 2D matrix of lagged variables. - */ +// Note that the return value is the value of the lagged order position, not the index. std::vector> CppLaggedVar4Grid( - std::vector> mat, - int lagNum) { - // Get the number of columns and rows in the input matrix - int numCols = mat[0].size(); + const std::vector>& mat, + int lagNum +) { int numRows = mat.size(); + int numCols = mat[0].size(); + int totalElements = numRows * numCols; + std::vector> result(totalElements, std::vector(8, std::numeric_limits::quiet_NaN())); - // Expand the matrix using the ExpandMat function - mat = ExpandMat(mat, lagNum); - - // Initialize the lagged variable matrix with dimensions (numRows * numCols) x (8 * lagNum) - int laggedVarRows = numRows * numCols; - int laggedVarCols = 8 * lagNum; - std::vector> laggedVar(laggedVarRows, std::vector(laggedVarCols, std::numeric_limits::quiet_NaN())); - - // Iterate over each row and column of the original matrix - for (int r = 1; r <= numRows; ++r) { - for (int c = 1; c <= numCols; ++c) { - int item = 0; // Index for the laggedVar matrix - - // Generate the sequence of column offsets - int colsAddStart = -std::trunc((3 + (lagNum - 1) * 2) / 2); - int colsAddEnd = std::trunc((3 + (lagNum - 1) * 2) / 2); - for (int la = colsAddStart; la <= colsAddEnd; ++la) { - laggedVar[(r - 1) * numCols + c - 1][item] = mat[r - 1 + lagNum][c - 1 + la + lagNum]; - item++; - } - - // Generate the sequence of row offsets - for (int ra = -(lagNum - 1); ra <= (lagNum - 1); ++ra) { - laggedVar[(r - 1) * numCols + c - 1][item] = mat[r - 1 + ra + lagNum][c - 1 - lagNum + lagNum]; - item++; - - laggedVar[(r - 1) * numCols + c - 1][item] = mat[r - 1 + ra + lagNum][c - 1 + lagNum + lagNum]; - item++; - } + // Directions for Moore neighborhood (8 directions) + const int dx[] = {-1, -1, -1, 0, 0, 1, 1, 1}; + const int dy[] = {-1, 0, 1, -1, 1, -1, 0, 1}; - // Generate the sequence of column offsets again - for (int la = colsAddStart; la <= colsAddEnd; ++la) { - laggedVar[(r - 1) * numCols + c - 1][item] = mat[r - 1 + lagNum + lagNum][c - 1 + la + lagNum]; - item++; + for (int r = 0; r < numRows; ++r) { + for (int c = 0; c < numCols; ++c) { + int elementIndex = r * numCols + c; + int neighborIndex = 0; + + // Calculate positions that are exactly 'lagNum' units away + for (int dir = 0; dir < 8; ++dir) { + int dr = dx[dir] * lagNum; + int dc = dy[dir] * lagNum; + int rNeighbor = r + dr; + int cNeighbor = c + dc; + + // Check if the neighbor position is within bounds + if (rNeighbor >= 0 && rNeighbor < numRows && cNeighbor >= 0 && cNeighbor < numCols) { + result[elementIndex][neighborIndex] = mat[rNeighbor][cNeighbor]; + } + // Else, it remains NaN + + neighborIndex++; + if (neighborIndex >= 8) { + break; + } } } } - return laggedVar; + return result; } -/** - * Generates a list of grid embeddings for the input matrix and embedding dimension. - * - * @param mat The input matrix represented as a 2D vector of doubles. - * @param E The embedding dimension, which determines the number of embeddings to generate. - * @return A vector of 2D vectors, where each element is an embedding of the input matrix. - */ std::vector>> GenGridEmbeddings( const std::vector>& mat, int E) { // Initialize a vector to store the embeddings diff --git a/src/CppGridUtils.h b/src/CppGridUtils.h index d2f7f46..52a74fc 100644 --- a/src/CppGridUtils.h +++ b/src/CppGridUtils.h @@ -7,11 +7,13 @@ #include std::vector> CppLaggedVar4Grid( - std::vector> mat, - int lagNum); + const std::vector>& mat, + int lagNum +); std::vector>> GenGridEmbeddings( const std::vector>& mat, - int E); + int E +); #endif // CppGridUtils_H diff --git a/src/CppLatticeUtils.cpp b/src/CppLatticeUtils.cpp index e7dae04..ca325dd 100644 --- a/src/CppLatticeUtils.cpp +++ b/src/CppLatticeUtils.cpp @@ -36,6 +36,7 @@ std::vector> nb2vec(Rcpp::List nb) { } // Function to compute lagged neighborhoods for a given lag number +// **Note that the return value corresponds to the cumulative neighbor indices for lagNum** std::vector> CppLaggedVar4Lattice(std::vector> spNeighbor, int lagNum) { // If lagNum is less than 1, return an empty vector diff --git a/src/GCCM4Grid.cpp b/src/GCCM4Grid.cpp new file mode 100644 index 0000000..fca5d0a --- /dev/null +++ b/src/GCCM4Grid.cpp @@ -0,0 +1,305 @@ +#include +#include +#include +#include +#include +#include +#include +#include "CppStats.h" +#include "CppGridUtils.h" +#include + +// Function to locate the index in a 2D grid +int locate(int curRow, int curCol, int totalRow, int totalCol) { + return (curRow - 1) * totalCol + curCol - 1; +} + +// Function to calculate the distance between embeddings, handling NaN values +std::vector DistCom(const std::vector>>& embeddings, + const std::vector& libs, int p) { + std::vector distances; + distances.reserve(libs.size()); + + // Iterate over each embedding + for (const auto& embedding : embeddings) { + // Get the p-th element of the current embedding + const std::vector& emd_p = embedding[p]; + + // Calculate the absolute difference between the p-th element and each library element + for (int lib : libs) { + const std::vector& emd_lib = embedding[lib]; + double distance = 0.0; + bool has_nan = false; // Flag to check if any NaN is encountered + + // Compute the distance while checking for NaN + for (size_t i = 0; i < emd_p.size(); ++i) { + if (std::isnan(emd_p[i]) || std::isnan(emd_lib[i])) { + has_nan = true; + break; + } + distance += std::abs(emd_p[i] - emd_lib[i]); + } + + // If NaN is encountered, skip this distance + if (has_nan) { + distances.push_back(std::numeric_limits::quiet_NaN()); + } else { + distances.push_back(distance); + } + } + } + + // Calculate row-wise mean, excluding NaN values + std::vector row_means(libs.size(), 0.0); + for (size_t i = 0; i < libs.size(); ++i) { + double sum = 0.0; + int count = 0; + + // Aggregate non-NaN distances + for (size_t e = 0; e < embeddings.size(); ++e) { + double distance = distances[i + e * libs.size()]; + if (!std::isnan(distance)) { + sum += distance; + ++count; + } + } + + // Compute mean, or set to NaN if all values are NaN + if (count > 0) { + row_means[i] = sum / count; + } else { + row_means[i] = std::numeric_limits::quiet_NaN(); + } + } + + return row_means; +} + +// Function to perform Simplex Projection Grid +double SimplexProjectionGrid(const std::vector>>& embeddings, + const std::vector& target, + const std::vector& lib_indices, + const std::vector& pred_indices, + int num_neighbors) { + std::vector pred(target.size(), std::numeric_limits::quiet_NaN()); + + for (size_t p = 0; p < pred_indices.size(); ++p) { + if (!pred_indices[p]) continue; + + // Temporarily exclude the current prediction index from the library + // bool temp_lib = lib_indices[p]; + std::vector modified_lib_indices = lib_indices; // Create a copy of lib_indices + modified_lib_indices[p] = false; + + // Get the indices of the library points + std::vector libs; + for (size_t i = 0; i < modified_lib_indices.size(); ++i) { + if (modified_lib_indices[i]) { + libs.push_back(i); + } + } + + // Calculate distances + std::vector distances = DistCom(embeddings, libs, p); + + // Find the nearest neighbors, excluding NaN distances + std::vector neighbors; + std::vector valid_distances; + for (size_t i = 0; i < distances.size(); ++i) { + if (!std::isnan(distances[i])) { + neighbors.push_back(i); + valid_distances.push_back(distances[i]); + } + } + + // If fewer than num_neighbors valid distances, skip this prediction + if (static_cast(neighbors.size()) < num_neighbors) { + continue; // Skip this prediction + } + + // Sort valid distances and get the top num_neighbors + std::vector sorted_indices(neighbors.size()); + std::iota(sorted_indices.begin(), sorted_indices.end(), 0); + std::sort(sorted_indices.begin(), sorted_indices.end(), [&](int a, int b) { + return valid_distances[a] < valid_distances[b]; + }); + + // Get the top num_neighbors + std::vector top_neighbors; + for (int i = 0; i < num_neighbors; ++i) { + top_neighbors.push_back(neighbors[sorted_indices[i]]); + } + + double min_distance = valid_distances[sorted_indices[0]]; + + // Compute weights + std::vector weights(num_neighbors, 0.000001); + if (min_distance == 0) { + for (int i = 0; i < num_neighbors; ++i) { + if (valid_distances[sorted_indices[i]] == 0) { + weights[i] = 1.0; + } + } + } else { + for (int i = 0; i < num_neighbors; ++i) { + weights[i] = std::exp(-valid_distances[sorted_indices[i]] / min_distance); + if (weights[i] < 0.000001) { + weights[i] = 0.000001; + } + } + } + + double total_weight = std::accumulate(weights.begin(), weights.end(), 0.0); + + // Make prediction + double prediction = 0.0; + for (int i = 0; i < num_neighbors; ++i) { + prediction += weights[i] * target[libs[top_neighbors[i]]]; + } + pred[p] = prediction / total_weight; + } + + // Return the Pearson correlation between the target and the predicted values + return PearsonCor(target, pred, true); +} + +// GCCMSingle4Grid function +std::vector> GCCMSingle4Grid( + const std::vector>>& xEmbedings, + const std::vector& yPred, + int lib_size, + const std::vector>& pred, + int totalRow, + int totalCol, + int b) { + + std::vector> x_xmap_y; + + for (int r = 1; r <= totalRow - lib_size + 1; ++r) { + for (int c = 1; c <= totalCol - lib_size + 1; ++c) { + + // Initialize prediction and library indices + std::vector pred_indices(totalRow * totalCol, false); + std::vector lib_indices(totalRow * totalCol, false); + + // Set prediction indices + for (const auto& p : pred) { + pred_indices[locate(p.first, p.second, totalRow, totalCol)] = true; + } + + // Exclude NA values in yPred from prediction indices + for (size_t i = 0; i < yPred.size(); ++i) { + if (std::isnan(yPred[i])) { + pred_indices[i] = false; + } + } + + // Set library indices + for (int i = r; i < r + lib_size; ++i) { + for (int j = c; j < c + lib_size; ++j) { + lib_indices[locate(i, j, totalRow, totalCol)] = true; + } + } + + // Check if more than half of the library is NA + int na_count = 0; + for (size_t i = 0; i < lib_indices.size(); ++i) { + if (lib_indices[i] && std::isnan(yPred[i])) { + ++na_count; + } + } + if (na_count > (lib_size * lib_size) / 2) { + continue; + } + + // Run cross map and store results + double results = SimplexProjectionGrid(xEmbedings, yPred, lib_indices, pred_indices, b); + x_xmap_y.emplace_back(lib_size, results); + } + } + + return x_xmap_y; +} + +// GCCM4Grid function +std::vector> GCCM4Grid( + const std::vector>& xMatrix, + const std::vector>& yMatrix, + const std::vector& lib_sizes, + const std::vector>& pred, + int E, + int tau = 1, + int b = 0) { + + // If b is not provided, default it to E + 1 + if (b == 0) { + b = E + 1; + } + + // Get the dimensions of the xMatrix + int totalRow = xMatrix.size(); + int totalCol = xMatrix[0].size(); + + // Flatten yMatrix into a 1D array (row-major order) + std::vector yPred; + for (const auto& row : yMatrix) { + yPred.insert(yPred.end(), row.begin(), row.end()); + } + + // Generate embeddings for xMatrix + auto xEmbedings = GenGridEmbeddings(xMatrix, E); + + // Ensure the minimum value in unique_lib_sizes is E + 2 and remove duplicates + std::vector unique_lib_sizes; + std::unique_copy(lib_sizes.begin(), lib_sizes.end(), std::back_inserter(unique_lib_sizes), + [&](int a, int b) { return a == b; }); + std::transform(unique_lib_sizes.begin(), unique_lib_sizes.end(), unique_lib_sizes.begin(), + [&](int size) { return std::max(size, E + 2); }); + + // Initialize the result container + std::vector> x_xmap_y; + + // Iterate over each library size + for (int lib_size : unique_lib_sizes) { + // Perform single grid cross-mapping for the current library size + auto results = GCCMSingle4Grid(xEmbedings, yPred, lib_size, pred, totalRow, totalCol, b); + + // Append the results to the main result container + x_xmap_y.insert(x_xmap_y.end(), results.begin(), results.end()); + } + + // Perform the operations using RcppThread + RcppThread::ProgressBar bar(unique_lib_sizes.size(), 1); + RcppThread::parallelFor(0, unique_lib_sizes.size(), [&](size_t i) { + int lib_size = unique_lib_sizes[i]; + auto results = GCCMSingle4Grid(xEmbedings, yPred, lib_size, pred, totalRow, totalCol, b); + x_xmap_y.insert(x_xmap_y.end(), results.begin(), results.end()); + bar++; + }); + + // Group by the first int and compute the mean + std::map> grouped_results; + for (const auto& result : x_xmap_y) { + grouped_results[result.first].push_back(result.second); + } + + std::vector> final_results; + for (const auto& group : grouped_results) { + double mean_value = CppMean(group.second, true); + final_results.push_back({static_cast(group.first), mean_value}); + } + + int n = pred.size(); + // Calculate significance and confidence interval for each result + for (size_t i = 0; i < final_results.size(); ++i) { + double rho = final_results[i][1]; + double significance = CppSignificance(rho, n); + std::vector confidence_interval = CppConfidence(rho, n); + + final_results[i].push_back(significance); + final_results[i].push_back(confidence_interval[0]); + final_results[i].push_back(confidence_interval[1]); + } + + return final_results; +} diff --git a/src/GCCM4Grid.h b/src/GCCM4Grid.h new file mode 100644 index 0000000..db0440f --- /dev/null +++ b/src/GCCM4Grid.h @@ -0,0 +1,42 @@ +#ifndef GCCM4Grid_H +#define GCCM4Grid_H + +#include +#include +#include +#include +#include +#include +#include +#include "CppStats.h" +#include "CppGridUtils.h" +#include + +// Function to perform Simplex Projection Grid +double SimplexProjectionGrid(const std::vector>>& embeddings, + const std::vector& target, + const std::vector& lib_indices, + const std::vector& pred_indices, + int num_neighbors); + +// GCCMSingle4Grid function +std::vector> GCCMSingle4Grid( + const std::vector>>& xEmbedings, + const std::vector& yPred, + int lib_size, + const std::vector>& pred, + int totalRow, + int totalCol, + int b); + +// GCCM4Grid function +std::vector> GCCM4Grid( + const std::vector>& xMatrix, + const std::vector>& yMatrix, + const std::vector& lib_sizes, + const std::vector>& pred, + int E, + int tau = 1, + int b = 0); + +#endif // GCCM4Grid_H diff --git a/src/GridExp.cpp b/src/GridExp.cpp index 599aba4..8148ab0 100644 --- a/src/GridExp.cpp +++ b/src/GridExp.cpp @@ -1,5 +1,6 @@ #include #include "CppGridUtils.h" +#include "GCCM4Grid.h" #include // [[Rcpp::export]] @@ -15,7 +16,7 @@ Rcpp::NumericMatrix RcppLaggedVar4Grid(Rcpp::NumericMatrix mat, int lagNum) { } } - // Call the Cpp function + // Call the CppLaggedVar4Grid function std::vector> laggedMat = CppLaggedVar4Grid(cppMat, lagNum); // Convert the result back to Rcpp::NumericMatrix @@ -67,3 +68,60 @@ Rcpp::List RcppGenGridEmbeddings(Rcpp::NumericMatrix mat, int E) { return result; } + +// [[Rcpp::export]] +Rcpp::NumericMatrix RcppGCCM4Grid( + const Rcpp::NumericMatrix& xMatrix, + const Rcpp::NumericMatrix& yMatrix, + const Rcpp::IntegerVector& lib_sizes, + const Rcpp::IntegerMatrix& pred, + int E) { + + // Convert Rcpp NumericMatrix to std::vector> + std::vector> xMatrix_cpp(xMatrix.nrow(), std::vector(xMatrix.ncol())); + for (int i = 0; i < xMatrix.nrow(); ++i) { + for (int j = 0; j < xMatrix.ncol(); ++j) { + xMatrix_cpp[i][j] = xMatrix(i, j); + } + } + + // Convert Rcpp NumericMatrix to std::vector> + std::vector> yMatrix_cpp(yMatrix.nrow(), std::vector(yMatrix.ncol())); + for (int i = 0; i < yMatrix.nrow(); ++i) { + for (int j = 0; j < yMatrix.ncol(); ++j) { + yMatrix_cpp[i][j] = yMatrix(i, j); + } + } + + // Convert Rcpp IntegerVector to std::vector + std::vector lib_sizes_cpp(lib_sizes.size()); + for (int i = 0; i < lib_sizes.size(); ++i) { + lib_sizes_cpp[i] = lib_sizes[i]; + } + + // Convert Rcpp IntegerMatrix to std::vector> + std::vector> pred_cpp(pred.nrow()); + for (int i = 0; i < pred.nrow(); ++i) { + pred_cpp[i] = std::make_pair(pred(i, 0), pred(i, 1)); + } + + // Call the C++ function GCCM4Grid + std::vector> result = GCCM4Grid( + xMatrix_cpp, + yMatrix_cpp, + lib_sizes_cpp, + pred_cpp, + E + ); + + Rcpp::NumericMatrix resultMatrix(result.size(), 5); + for (size_t i = 0; i < result.size(); ++i) { + resultMatrix(i, 0) = result[i][0]; + resultMatrix(i, 1) = result[i][1]; + resultMatrix(i, 2) = result[i][2]; + resultMatrix(i, 3) = result[i][3]; + resultMatrix(i, 4) = result[i][4]; + } + + return resultMatrix; +} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 6b30431..39fdda3 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -35,6 +35,21 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// RcppGCCM4Grid +Rcpp::NumericMatrix RcppGCCM4Grid(const Rcpp::NumericMatrix& xMatrix, const Rcpp::NumericMatrix& yMatrix, const Rcpp::IntegerVector& lib_sizes, const Rcpp::IntegerMatrix& pred, int E); +RcppExport SEXP _spEDM_RcppGCCM4Grid(SEXP xMatrixSEXP, SEXP yMatrixSEXP, SEXP lib_sizesSEXP, SEXP predSEXP, SEXP ESEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const Rcpp::NumericMatrix& >::type xMatrix(xMatrixSEXP); + Rcpp::traits::input_parameter< const Rcpp::NumericMatrix& >::type yMatrix(yMatrixSEXP); + Rcpp::traits::input_parameter< const Rcpp::IntegerVector& >::type lib_sizes(lib_sizesSEXP); + Rcpp::traits::input_parameter< const Rcpp::IntegerMatrix& >::type pred(predSEXP); + Rcpp::traits::input_parameter< int >::type E(ESEXP); + rcpp_result_gen = Rcpp::wrap(RcppGCCM4Grid(xMatrix, yMatrix, lib_sizes, pred, E)); + return rcpp_result_gen; +END_RCPP +} // RcppConfidence Rcpp::NumericVector RcppConfidence(double r, int n, double level); RcppExport SEXP _spEDM_RcppConfidence(SEXP rSEXP, SEXP nSEXP, SEXP levelSEXP) { @@ -106,6 +121,7 @@ END_RCPP static const R_CallMethodDef CallEntries[] = { {"_spEDM_RcppLaggedVar4Grid", (DL_FUNC) &_spEDM_RcppLaggedVar4Grid, 2}, {"_spEDM_RcppGenGridEmbeddings", (DL_FUNC) &_spEDM_RcppGenGridEmbeddings, 2}, + {"_spEDM_RcppGCCM4Grid", (DL_FUNC) &_spEDM_RcppGCCM4Grid, 5}, {"_spEDM_RcppConfidence", (DL_FUNC) &_spEDM_RcppConfidence, 3}, {"_spEDM_RcppLinearTrendRM", (DL_FUNC) &_spEDM_RcppLinearTrendRM, 4}, {"_spEDM_RcppLaggedVar4Lattice", (DL_FUNC) &_spEDM_RcppLaggedVar4Lattice, 2}, diff --git a/vignettes/GCCM.Rmd.orig b/vignettes/GCCM.Rmd.orig index 622a18a..5234a5f 100644 --- a/vignettes/GCCM.Rmd.orig +++ b/vignettes/GCCM.Rmd.orig @@ -68,12 +68,16 @@ $$ z = \frac{1}{2} \ln \left(\frac{1+\rho}{1-\rho}\right) $$ -## 2. An example of county-level population density +## 2. Examples + +### 2.1 Install the `spEDM` package ```r install.packages("spEDM", dep = TRUE) ``` +### 2.2 An example of lattice data about county-level population density + Load data and package: ```{r} @@ -92,25 +96,25 @@ Run GCCM: ```{r} startTime = Sys.time() -res = gccm(cause = "Pre", - effect = "popDensity", - nb = popd_nb, - coords = c("x","y"), - data = popdensity, - libsizes = seq(10, 2800, by = 100), - E = 3) +pd_res = gccm(cause = "Pre", + effect = "popDensity", + nb = popd_nb, + coords = c("x","y"), + data = popdensity, + libsizes = seq(10, 2800, by = 100), + E = 3) endTime = Sys.time() print(difftime(endTime,startTime, units ="mins")) -res +pd_res ``` Visualize the result: -```{r gccm_plot, fig.width=6.75,fig.height=5.5,fig.dpi=300,fig.cap=knitr::asis_output("**Figure 1**. The cross-mapping prediction outputs between population density and county-level Precipitation.")} +```{r fig1,fig.width=6.75,fig.height=5.5,fig.dpi=300,fig.cap=knitr::asis_output("**Figure 1**. The cross-mapping prediction outputs between population density and county-level Precipitation.")} windowsFonts(TNR = windowsFont("Times New Roman")) -fig_gccm = ggplot2::ggplot(data = res, - ggplot2::aes(x = lib_sizes)) + +fig1 = ggplot2::ggplot(data = pd_res, + ggplot2::aes(x = lib_sizes)) + ggplot2::geom_line(ggplot2::aes(y = x_xmap_y_mean, color = "x xmap y"), lwd = 1.25) + @@ -137,6 +141,37 @@ fig_gccm = ggplot2::ggplot(data = res, legend.justification = c('right','top'), legend.background = ggplot2::element_rect(fill = 'transparent'), legend.text = ggplot2::element_text(family = "TNR")) -fig_gccm +fig1 ``` +### 2.3 An example of grid data about soil pollution + +Load data and package: + +```{r fig2,fig.width=6.55,fig.height=2.60,fig.dpi=300,fig.cap=knitr::asis_output("**Figure 2**. Maps of soil Cu concentration, the density of industrial pollutants and nightlight.")} +library(spEDM) + +soilpollution = terra::rast(system.file("extdata/soilpollution.tif", + package = "spEDM")) +soilpollution + +terra::plot(soilpollution,nc = 3, + mar = rep(0.1,4), + oma = rep(0.1,4), + axes = FALSE, + legend = FALSE) +``` + +Run GCCM: + +```{r} +startTime = Sys.time() +sp_res = gccm(cause = "industry", + effect = "cu", + data = soilpollution, + libsizes = seq(10,120,20), + E = 3) +endTime = Sys.time() +print(difftime(endTime,startTime, units ="mins")) +sp_res +```