draft metacoupling support

NAMESPACE

@@ -1,4 +1,5 @@
 # Generated by roxygen2: do not edit by hand
 
 export(ccd)
+export(metacoupling)
 useDynLib(coupling, .registration = TRUE)

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

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

_pkgdown.yml

@@ -30,3 +30,7 @@ reference:
 - subtitle: coupling coordination degree
   contents:
   - ccd
+
+- subtitle: meta-coupling analysis
+  contents:
+  - metacoupling

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)

inst/include/coupling/ccd.hpp

@@ -178,6 +178,7 @@
 #include <algorithm>
 #include <stdexcept>
 #include <RcppThread.h>
+#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<double>::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());