diff --git a/NAMESPACE b/NAMESPACE index 7f5f347..bc23093 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,4 +1,5 @@ # Generated by roxygen2: do not edit by hand export(ccd) +export(metacoupling) useDynLib(coupling, .registration = TRUE) diff --git a/R/RcppExports.R b/R/RcppExports.R index bcf759d..5f41588 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -5,3 +5,7 @@ RcppCCD <- function(mat, weight, method = "standard", threads = 1L) { .Call(`_coupling_RcppCCD`, mat, weight, method, threads) } +RcppMetaCoupling <- function(mat, swm_peri, swm_tele, weight, method = "standard", threads = 1L) { + .Call(`_coupling_RcppMetaCoupling`, mat, swm_peri, swm_tele, weight, method, threads) +} + diff --git a/R/metacoupling.R b/R/metacoupling.R new file mode 100644 index 0000000..9168d51 --- /dev/null +++ b/R/metacoupling.R @@ -0,0 +1,44 @@ +#' Metacoupling Analysis +#' +#' @inheritParams ccd +#' @param swm_peri A numeric matrix representing the **peri (local) spatial weight matrix**. +#' Must be square with dimension equal to `nrow(data)`. If `NULL`, a zero matrix is used. +#' @param swm_tele A numeric matrix representing the **tele (long-distance) spatial weight matrix**. +#' Must be square with dimension equal to `nrow(data)`. If `NULL`, a zero matrix is used. +#' +#' @return A data.frame with: +#' \itemize{ +#' \item `Intra_C`: intra-system coupling degree +#' \item `Intra_D`: intra-system coordination degree +#' \item `Peri_C`: peri-coupling degree +#' \item `Peri_D`: peri coordination degree +#' \item `Tele_C`: tele-coupling degree +#' \item `Tele_D`: tele coordination degree +#' } +#' +#' @export +#' +#' @details +#' Full model definitions and formulas are available at: +#' \url{https://github.com/stscl/coupling/discussions/8} +#' +#' @note +#' Input values should be normalized to `[0, 1]`. Spatial weight matrices are +#' typically symmetric. +#' +#' @examples +#' set.seed(42) +#' mat = matrix(runif(20), nrow = 5) +#' swm1 = apply(matrix(runif(25), 5, 5), 1, \(.x) .x / sum(.x)) +#' swm2 = apply(matrix(runif(25), 5, 5), 1, \(.x) .x / sum(.x)) +#' coupling::metacoupling(mat, swm1, swm2) +#' +metacoupling = \(data, swm_peri = NULL, swm_tele = NULL, weight = NULL, + method = c("standard", "wang", "fan"), threads = 1){ + mat = as.matrix(data) + method = match.arg(method) + if (is.null(weight)) weight = rep(1, times = ncol(mat)) / ncol(mat) + if (is.null(swm_peri)) swm_peri = matrix(0, nrow(data), nrow(data)) + if (is.null(swm_tele)) swm_tele = matrix(0, nrow(data), nrow(data)) + return(RcppMetaCoupling(mat, swm_peri, swm_tele, weight, method, threads)) +} diff --git a/_pkgdown.yml b/_pkgdown.yml index b3d989a..36d1273 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -30,3 +30,7 @@ reference: - subtitle: coupling coordination degree contents: - ccd + +- subtitle: meta-coupling analysis + contents: + - metacoupling diff --git a/inst/include/coupling.h b/inst/include/coupling.h index 0456ec5..6d73838 100644 --- a/inst/include/coupling.h +++ b/inst/include/coupling.h @@ -39,6 +39,7 @@ #include "coupling/numericutils.hpp" #include "coupling/ccd.hpp" +#include "coupling/metacoupling.hpp" // ============================================================ // Convenience Converters (Inline helpers for R/C++ interop) diff --git a/inst/include/coupling/ccd.hpp b/inst/include/coupling/ccd.hpp index e605273..d978cb8 100644 --- a/inst/include/coupling/ccd.hpp +++ b/inst/include/coupling/ccd.hpp @@ -178,6 +178,7 @@ #include #include #include +#include "coupling/numericutils.hpp" namespace coupling { @@ -196,6 +197,7 @@ inline double ccd_c_single( const std::string& method = "standard" ) { size_t p = vec.size(); // number of U values + if (p <= 1) return std::numeric_limits::quiet_NaN(); double C_val = 0.0; @@ -233,7 +235,7 @@ inline double ccd_c_single( double denom = (p - 1) * p / 2.0; double term1 = 1.0 - (sum_dist / denom); - if (term1 < 0) term1 = 0; + // if (term1 < 0) term1 = 0; double max_u = *std::max_element(vec.begin(), vec.end()); diff --git a/inst/include/coupling/metacoupling.hpp b/inst/include/coupling/metacoupling.hpp new file mode 100644 index 0000000..5e0b639 --- /dev/null +++ b/inst/include/coupling/metacoupling.hpp @@ -0,0 +1,454 @@ +/****************************************************************************** + * File: metacoupling.hpp + * + * Metacoupling Model (Extension of CCD) + * ------------------------------------------------ + * + * This module implements the metacoupling framework, an extension of the + * Coupling Coordination Degree (CCD) model, designed to measure both + * intra-system and inter-system coupling relationships across spatial units. + * + * The metacoupling model integrates: + * + * • intra-unit coupling (within a spatial unit) + * • peri-coupling (short-distance interactions) + * • tele-coupling (long-distance interactions) + * + * It further extends the CCD framework by introducing a permutation-based + * cross-system coupling mechanism. + * + * --------------------------------------------------------------------------- + * Mathematical formulation + * --------------------------------------------------------------------------- + * + * (See https://github.com/stscl/coupling/discussions/8 for a better reading experience.) + * + * Let there be: + * + * • n spatial units + * • p indicators per unit + * + * For unit i: + * + * U_i = {U_{i1}, U_{i2}, ..., U_{ip}} + * + * --------------------------------------------------------------------------- + * 1. Intra-system coupling (C_intra) + * --------------------------------------------------------------------------- + * + * Standard CCD coupling computed for each unit: + * + * C_intra(i) = CCD(U_i) + * + * where CCD is defined in ccd.hpp. + * + * --------------------------------------------------------------------------- + * 2. Inter-system metacoupling (C_ij) + * --------------------------------------------------------------------------- + * + * For two units i and j, metacoupling is defined by constructing + * mixed systems through all possible indicator-level combinations: + * + * V(mask) = {v_1, v_2, ..., v_p} + * + * where: + * + * v_k = U_{ik} if bit k in mask = 1 + * v_k = U_{jk} if bit k in mask = 0 + * + * All binary masks are enumerated: + * + * mask ∈ {0, 1}^p + * + * EXCLUDING: + * + * • mask = 0 (all indicators from j) + * • mask = 2^p − 1 (all indicators from i) + * + * Total combinations used: + * + * 2^p − 2 + * + * The inter-system coupling is: + * + * C_ij = mean over all valid masks of CCD(V(mask)) + * + * --------------------------------------------------------------------------- + * 3. Spatial aggregation (peri / tele) + * --------------------------------------------------------------------------- + * + * Given spatial weight matrices: + * + * W_peri (short-distance interactions) + * W_tele (long-distance interactions) + * + * Both are assumed to be symmetric. + * + * Aggregated coupling: + * + * C_peri(i) = Σ_j W_peri(i,j) * C_ij + * C_tele(i) = Σ_j W_tele(i,j) * C_ij + * + * --------------------------------------------------------------------------- + * 4. Composite development index (T) + * --------------------------------------------------------------------------- + * + * For each system: + * + * T = Σ_k w_k * U_k + * + * For inter-system: + * + * T_ij = mean over all permutations of T(V(mask)) + * + * --------------------------------------------------------------------------- + * 5. Coupling coordination degree (D) + * --------------------------------------------------------------------------- + * + * Final metric: + * + * D = sqrt(C × T) + * + * Applied to: + * + * • intra-system + * • peri-coupling + * • tele-coupling + * + * --------------------------------------------------------------------------- + * Input data format + * --------------------------------------------------------------------------- + * + * mat : vector> + * + * • rows = spatial units + * • columns = indicators + * • values MUST be normalized to [0,1] + * + * swm_peri / swm_tele: + * + * • symmetric spatial weight matrices (n × n) + * + * weight: + * + * • length = number of indicators + * + * --------------------------------------------------------------------------- + * Output + * --------------------------------------------------------------------------- + * + * metacoupling_c: + * + * result is a 3 × n_units matrix: + * + * result[0][i] = C_intra(i) + * result[1][i] = C_peri(i) + * result[2][i] = C_tele(i) + * + * metacoupling: + * + * result is a 6 × n_units matrix: + * + * result[0][i] = C_intra(i) + * result[1][i] = D_intra(i) + * result[2][i] = C_peri(i) + * result[3][i] = D_peri(i) + * result[4][i] = C_tele(i) + * result[5][i] = D_tele(i) + * + * --------------------------------------------------------------------------- + * Notes + * --------------------------------------------------------------------------- + * + * • The number of inter-system combinations grows exponentially (2^p − 2) + * + * • This framework captures cross-unit structural interactions + * at the indicator level, rather than treating systems as indivisible units + * + * • Input normalization is critical for meaningful interpretation + * + * • Computational cost increases rapidly with p + * + * --------------------------------------------------------------------------- + * Author: Wenbo Lyu (Github: @SpatLyu) + * License: GPL-3 + ******************************************************************************/ + +#ifndef COUPLING_METACOUPLING_HPP +#define COUPLING_METACOUPLING_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include "coupling/ccd.hpp" +#include "coupling/numericutils.hpp" + +namespace coupling +{ + +namespace metacoupling +{ + +inline std::vector> metacoupling_c( + const std::vector>& mat, + const std::vector>& swm_peri, + const std::vector>& swm_tele, + const std::string& method = "standard", + size_t threads = 1 +) { + // Number of spatial units (observations) + size_t n_units = mat.size(); + + // Number of indicators (system dimensions) + size_t p = mat[0].size(); + + // Output: + // [0][i] = intra-coupling (original CCD for unit i) + // [1][i] = peri-coupling + // [2][i] = tele-coupling + std::vector> result( + 3, + std::vector(n_units, std::numeric_limits::quiet_NaN()) + ); + + if (p <= 1) return result; + + // Total number of binary combinations = 2^p + size_t full_perm = 1ULL << p; + + // ============================================================ + // Compute metacoupling between unit i and j + // + // Definition: + // - Enumerate ALL mixed systems + // - Each dimension independently comes from i or j + // - EXCLUDE: + // mask == 0 → all from j + // mask == full_perm-1 → all from i + // + // Total combinations used = 2^p - 2 + // ============================================================ + auto compute_pair = [&](size_t i, size_t j) { + + double sum_c = 0.0; + size_t count = 0; + + for (size_t mask = 0; mask < full_perm; ++mask) { + + // ---------------------------------------------------- + // Remove trivial (non-mixed) systems + // ---------------------------------------------------- + if (mask == 0 || mask == (full_perm - 1)) continue; + + // ---------------------------------------------------- + // Construct mixed system vector + // + // Bit k of mask determines source: + // 1 → take from unit i + // 0 → take from unit j + // + // Example (p = 3): + // mask = 101 → [i, j, i] + // ---------------------------------------------------- + std::vector vec(p); + + for (size_t k = 0; k < p; ++k) { + if (mask & (1ULL << k)) { + vec[k] = mat[i][k]; + } else { + vec[k] = mat[j][k]; + } + } + + // Compute CCD coupling + sum_c += coupling::ccd::ccd_c_single(vec, method); + ++count; + } + + // count should be exactly 2^p - 2 + return sum_c / static_cast(count); + }; + + // ============================================================ + // Main worker for each spatial unit s + // ============================================================ + auto worker = [&](size_t s) { + + // 1. Intra-coupling (original system) + result[0][s] = coupling::ccd::ccd_c_single(mat[s], method); + + double peri_sum = 0.0; + double tele_sum = 0.0; + + // 2. Iterate over all other units + for (size_t j = 0; j < n_units; ++j) { + + double w_peri = swm_peri[s][j]; + double w_tele = swm_tele[s][j]; + + // Skip if no interaction + if (coupling::numericutils::doubleNearlyEqual(w_peri, 0.0) && + coupling::numericutils::doubleNearlyEqual(w_tele, 0.0)) continue; + + // Compute symmetric coupling between (s, j) + double c_val = compute_pair(s, j); + + // 3. Weighted aggregation + // Since SWM is symmetric, no need to compute (j, s) separately + if (!coupling::numericutils::doubleNearlyEqual(w_peri, 0.0)) { + peri_sum += w_peri * c_val; + } + + if (!coupling::numericutils::doubleNearlyEqual(w_tele, 0.0)) { + tele_sum += w_tele * c_val; + } + } + + result[1][s] = peri_sum; + result[2][s] = tele_sum; + }; + + // ============================================================ + // Parallel or serial execution + // ============================================================ + if (threads <= 1) { + for (size_t s = 0; s < n_units; ++s) { + worker(s); + } + } else { + RcppThread::parallelFor(0, n_units, worker, threads); + } + + return result; +} + +inline std::vector> metacoupling( + const std::vector>& mat, + const std::vector>& swm_peri, + const std::vector>& swm_tele, + const std::vector& weight, + const std::string& method = "standard", + size_t threads = 1 +) { + size_t n_units = mat.size(); + size_t p = mat[0].size(); + + // ============================================================ + // Output: 6 × n_units + // + // [0] = Intra_C + // [1] = Intra_D + // [2] = Peri_C + // [3] = Peri_D + // [4] = Tele_C + // [5] = Tele_D + // ============================================================ + std::vector> result( + 6, + std::vector(n_units, std::numeric_limits::quiet_NaN()) + ); + + if (p <= 1) return result; + + size_t full_perm = 1ULL << p; + + // ============================================================ + // Step 1: reuse existing C computation + // ============================================================ + auto C_res = metacoupling_c(mat, swm_peri, swm_tele, method, threads); + + // ============================================================ + // Step 2: T computation + // ============================================================ + auto compute_T_intra = [&](size_t i) { + double t = 0.0; + for (size_t k = 0; k < p; ++k) { + t += weight[k] * mat[i][k]; + } + return t; + }; + + auto compute_pair_T = [&](size_t i, size_t j) { + + double sum_t = 0.0; + size_t count = 0; + + for (size_t mask = 0; mask < full_perm; ++mask) { + + if (mask == 0 || mask == (full_perm - 1)) continue; + + double t_val = 0.0; + + for (size_t k = 0; k < p; ++k) { + if (mask & (1ULL << k)) { + t_val += weight[k] * mat[i][k]; + } else { + t_val += weight[k] * mat[j][k]; + } + } + + sum_t += t_val; + ++count; + } + + return sum_t / static_cast(count); + }; + + // ============================================================ + // Step 3: assemble results + // ============================================================ + for (size_t s = 0; s < n_units; ++s) { + + // ---- intra ---- + double C_intra = C_res[0][s]; + double T_intra = compute_T_intra(s); + + result[0][s] = C_intra; + result[1][s] = std::sqrt(C_intra * T_intra); + + double peri_T_sum = 0.0; + double tele_T_sum = 0.0; + + // ---- interactions ---- + for (size_t j = 0; j < n_units; ++j) { + + double w_peri = swm_peri[s][j]; + double w_tele = swm_tele[s][j]; + + if (coupling::numericutils::doubleNearlyEqual(w_peri, 0.0) && + coupling::numericutils::doubleNearlyEqual(w_tele, 0.0)) continue; + + double t_val = compute_pair_T(s, j); + + if (!coupling::numericutils::doubleNearlyEqual(w_peri, 0.0)) { + peri_T_sum += w_peri * t_val; + } + + if (!coupling::numericutils::doubleNearlyEqual(w_tele, 0.0)) { + tele_T_sum += w_tele * t_val; + } + } + + double C_peri = C_res[1][s]; + double C_tele = C_res[2][s]; + + result[2][s] = C_peri; + result[3][s] = std::sqrt(C_peri * peri_T_sum); + + result[4][s] = C_tele; + result[5][s] = std::sqrt(C_tele * tele_T_sum); + } + + return result; +} + +} // namespace metacoupling + +} + +#endif // COUPLING_METACOUPLING_HPP diff --git a/man/metacoupling.Rd b/man/metacoupling.Rd new file mode 100644 index 0000000..f318dee --- /dev/null +++ b/man/metacoupling.Rd @@ -0,0 +1,62 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/metacoupling.R +\name{metacoupling} +\alias{metacoupling} +\title{Metacoupling Analysis} +\usage{ +metacoupling( + data, + swm_peri = NULL, + swm_tele = NULL, + weight = NULL, + method = c("standard", "wang", "fan"), + threads = 1 +) +} +\arguments{ +\item{data}{A numeric matrix or data.frame. Rows are observations, +columns are indicators.} + +\item{swm_peri}{A numeric matrix representing the \strong{peri (local) spatial weight matrix}. +Must be square with dimension equal to \code{nrow(data)}. If \code{NULL}, a zero matrix is used.} + +\item{swm_tele}{A numeric matrix representing the \strong{tele (long-distance) spatial weight matrix}. +Must be square with dimension equal to \code{nrow(data)}. If \code{NULL}, a zero matrix is used.} + +\item{weight}{Numeric vector of indicator weights. Must have length equal +to \code{ncol(data)}. If \code{NULL}, equal weights are used.} + +\item{method}{Coupling model. One of \code{"standard"}, \code{"wang"}, or \code{"fan"}.} + +\item{threads}{Number of threads used in computation.} +} +\value{ +A data.frame with: +\itemize{ +\item \code{Intra_C}: intra-system coupling degree +\item \code{Intra_D}: intra-system coordination degree +\item \code{Peri_C}: peri-coupling degree +\item \code{Peri_D}: peri coordination degree +\item \code{Tele_C}: tele-coupling degree +\item \code{Tele_D}: tele coordination degree +} +} +\description{ +Metacoupling Analysis +} +\details{ +Full model definitions and formulas are available at: +\url{https://github.com/stscl/coupling/discussions/8} +} +\note{ +Input values should be normalized to \verb{[0, 1]}. Spatial weight matrices are +typically symmetric. +} +\examples{ +set.seed(42) +mat = matrix(runif(20), nrow = 5) +swm1 = apply(matrix(runif(25), 5, 5), 1, \(.x) .x / sum(.x)) +swm2 = apply(matrix(runif(25), 5, 5), 1, \(.x) .x / sum(.x)) +coupling::metacoupling(mat, swm1, swm2) + +} diff --git a/src/CCD.cpp b/src/CCD.cpp deleted file mode 100644 index 309c353..0000000 --- a/src/CCD.cpp +++ /dev/null @@ -1,27 +0,0 @@ -#include -#include -#include -#include -#include -#include "coupling.h" - -// Wrapper function to preform coupling coordination degree analysis -// [[Rcpp::export(rng = false)]] -Rcpp::DataFrame RcppCCD(const Rcpp::NumericMatrix& mat, - const Rcpp::NumericVector& weight, - const std::string& method = "standard", - int threads = 1) -{ - std::vector> mat_std = coupling::convert::mat_r2std(mat, true); - std::vector wt_std = Rcpp::as>(weight); - - // Call ccd function - std::vector> res = coupling::ccd::ccd( - mat_std, wt_std, method, static_cast(std::abs(threads)) - ); - - return Rcpp::DataFrame::create( - Rcpp::Named("C") = res[0], - Rcpp::Named("D") = res[1] - ); -} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index afb6ae8..8dae5df 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -25,9 +25,25 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// RcppMetaCoupling +Rcpp::DataFrame RcppMetaCoupling(const Rcpp::NumericMatrix& mat, const Rcpp::NumericMatrix& swm_peri, const Rcpp::NumericMatrix& swm_tele, const Rcpp::NumericVector& weight, const std::string& method, int threads); +RcppExport SEXP _coupling_RcppMetaCoupling(SEXP matSEXP, SEXP swm_periSEXP, SEXP swm_teleSEXP, SEXP weightSEXP, SEXP methodSEXP, SEXP threadsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::traits::input_parameter< const Rcpp::NumericMatrix& >::type mat(matSEXP); + Rcpp::traits::input_parameter< const Rcpp::NumericMatrix& >::type swm_peri(swm_periSEXP); + Rcpp::traits::input_parameter< const Rcpp::NumericMatrix& >::type swm_tele(swm_teleSEXP); + Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type weight(weightSEXP); + Rcpp::traits::input_parameter< const std::string& >::type method(methodSEXP); + Rcpp::traits::input_parameter< int >::type threads(threadsSEXP); + rcpp_result_gen = Rcpp::wrap(RcppMetaCoupling(mat, swm_peri, swm_tele, weight, method, threads)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_coupling_RcppCCD", (DL_FUNC) &_coupling_RcppCCD, 4}, + {"_coupling_RcppMetaCoupling", (DL_FUNC) &_coupling_RcppMetaCoupling, 6}, {NULL, NULL, 0} }; diff --git a/src/wrappers.cpp b/src/wrappers.cpp new file mode 100644 index 0000000..6173dcd --- /dev/null +++ b/src/wrappers.cpp @@ -0,0 +1,57 @@ +#include +#include +#include +#include +#include +#include "coupling.h" + +// Wrapper function to preform coupling coordination degree analysis +// [[Rcpp::export(rng = false)]] +Rcpp::DataFrame RcppCCD(const Rcpp::NumericMatrix& mat, + const Rcpp::NumericVector& weight, + const std::string& method = "standard", + int threads = 1) +{ + std::vector> mat_std = coupling::convert::mat_r2std(mat, true); + std::vector wt_std = Rcpp::as>(weight); + + // Call ccd function + std::vector> res = coupling::ccd::ccd( + mat_std, wt_std, method, static_cast(std::abs(threads)) + ); + + return Rcpp::DataFrame::create( + Rcpp::Named("C") = res[0], + Rcpp::Named("D") = res[1] + ); +} + +// Wrapper function to preform meta-coupling analysis +// [[Rcpp::export(rng = false)]] +Rcpp::DataFrame RcppMetaCoupling( + const Rcpp::NumericMatrix& mat, + const Rcpp::NumericMatrix& swm_peri, + const Rcpp::NumericMatrix& swm_tele, + const Rcpp::NumericVector& weight, + const std::string& method = "standard", + int threads = 1) +{ + std::vector> mat_std = coupling::convert::mat_r2std(mat, true); + std::vector> swm_peri_std = coupling::convert::mat_r2std(swm_peri, true); + std::vector> swm_tele_std = coupling::convert::mat_r2std(swm_tele, true); + std::vector wt_std = Rcpp::as>(weight); + + // Call metacoupling function + std::vector> res = coupling::metacoupling::metacoupling( + mat_std, swm_peri_std, swm_tele_std, wt_std, method, static_cast(std::abs(threads)) + ); + + return Rcpp::DataFrame::create( + Rcpp::Named("Intra_C") = res[0], + Rcpp::Named("Intra_D") = res[1], + Rcpp::Named("Peri_C") = res[2], + Rcpp::Named("Peri_D") = res[3], + Rcpp::Named("Tele_C") = res[4], + Rcpp::Named("Tele_D") = res[5] + ); +}