R/internal_dcor_helpers.R

In dcorBSS: Distance-Correlation Based Methods for Blind Source Separation and Dependence Analysis

# Internal helpers for dcorBSS

validate_matrix_input <- function(X) {
  X <- as.matrix(X)
  if (!is.numeric(X)) {
    stop("`X` must be numeric.", call. = FALSE)
  }
  if (nrow(X) < 2L) {
    stop("`X` must contain at least two observations.", call. = FALSE)
  }
  if (ncol(X) < 2L) {
    stop("`X` must contain at least two variables.", call. = FALSE)
  }
  X
}

my_whitener <- function(X, center, scatter) {
  X <- as.matrix(X)
  center <- as.numeric(center)
  scatter <- as.matrix(scatter)
  
  if (length(center) != ncol(X)) {
    stop("`center` must have length ncol(X).", call. = FALSE)
  }
  if (!all(dim(scatter) == c(ncol(X), ncol(X)))) {
    stop("`scatter` must be a square matrix with dimension ncol(X).",
         call. = FALSE)
  }
  if (!isTRUE(all.equal(scatter, t(scatter), check.attributes = FALSE))) {
    stop("`scatter` must be symmetric.", call. = FALSE)
  }
  
  eig <- eigen(scatter, symmetric = TRUE)
  if (any(eig$values <= 0)) {
    stop("`scatter` must be positive definite for whitening.", call. = FALSE)
  }
  
  whitening_matrix <- eig$vectors %*%
    diag(1 / sqrt(eig$values), nrow = ncol(X))
  X_centered <- sweep(X, 2L, center, FUN = "-")
  Z <- X_centered %*% whitening_matrix
  K <- t(whitening_matrix)
  
  list(K = K, Z = Z)
}

givens_rotation <- function(angle, p, index) {
  Q <- diag(p)
  cval <- cos(angle)
  sval <- sin(angle)
  
  i <- index[1L]
  j <- index[2L]
  
  Q[i, i] <- cval
  Q[j, j] <- cval
  Q[i, j] <- -sval
  Q[j, i] <- sval
  Q
}

theta_to_W <- function(theta) {
  p <- (1 + sqrt(1 + 8 * length(theta))) / 2
  if (p != floor(p)) {
    stop("`theta` must have length p(p - 1) / 2 for some integer p.",
         call. = FALSE)
  }
  
  p <- as.integer(p)
  W <- diag(p)
  idx <- 1L
  
  for (i in seq_len(p - 1L)) {
    for (j in (i + 1L):p) {
      Q <- givens_rotation(theta[idx], p = p, index = c(j, i))
      W <- Q %*% W
      idx <- idx + 1L
    }
  }
  
  W
}

current_objective_sequential_dcorICA <- function(thetas,
                                                 Z,
                                                 transform,
                                                 alpha = 0.9975,
                                                 block_size = NULL) {
  Z <- as.matrix(Z)
  
  if ((ncol(Z) - 1L) != length(thetas)) {
    stop("Wrong dimensions in `current_objective_sequential_dcorICA()`. ",
         call. = FALSE)
  }
  
  d <- ncol(Z)
  Q <- diag(d)
  idx <- 2L
  
  for (theta in thetas) {
    Q <- givens_rotation(theta, d, c(idx, 1L)) %*% Q
    idx <- idx + 1L
  }
  
  S <- tcrossprod(Z, Q)
  
  if (identical(transform, "none")) {
    x_part <- S[, 1L, drop = FALSE]
    y_part <- S[, -1L, drop = FALSE]
    value <- dcor_large(x_part, y_part, block_size = block_size)^2
  } else if (identical(transform, "bowl")) {
    x_part <- bowl_transform(
      matrix(S[, 1L], ncol = 1L),
      alpha = alpha,
      do_scale = FALSE
    )
    y_part <- bowl_transform(
      S[, -1L, drop = FALSE],
      alpha = alpha,
      do_scale = FALSE
    )
    value <- dcor_large(x_part, y_part, block_size = block_size)^2
  } else {
    stop("`transform` must be one of 'none' or 'bowl'.", call. = FALSE)
  }
  
  value
}

optimize_sequential_dcorICA <- function(Z,
                                        transform = "none",
                                        alpha = 0.9975,
                                        sweeps = 1L,
                                        seed = NULL,
                                        block_size = NULL,
                                        optimizer = "cobyla",
                                        control = list()) {
  if (!is.null(seed)) {
    set.seed(seed)
  }
  
  if (is.null(control)) control <- list()
  
  if (optimizer == "cobyla"){
    default_control <- list()
  } else if (optimizer == "bobyqa"){
    default_control <- list(
      maxfun = 10000,
      rhobeg = 0.2,
      rhoend = 1e-6,
      iprint = 0
    )
  } else {
    stop("Invalid optimizer chosen.")
  }
  control <- utils::modifyList(default_control, control)
  
  p <- ncol(Z)
  W_total <- diag(p)
  convergence_log <- vector("list", length = sweeps)
  
  for (sweep_idx in seq_len(sweeps)) {
    perm <- sample.int(p)
    inverse_perm <- order(perm)
    Z_current <- Z[, perm, drop = FALSE]
    
    current_index <- 1L
    thetas <- numeric(p * (p - 1L) / 2L)
    sweep_log <- vector("list", length = p - 1L)
    
    for (i in seq_len(p - 1L)) {
      n_theta <- p - i
      
      if (optimizer == "cobyla"){
        res <- nloptr::cobyla(
          x0 = rep(0, n_theta),
          fn = current_objective_sequential_dcorICA,
          lower = rep(-pi, n_theta),
          upper = rep(pi, n_theta),
          Z = Z_current,
          control = control,
          transform = transform,
          alpha = alpha,
          block_size = block_size
        )
      } else if (optimizer == "bobyqa"){
        res <- minqa::bobyqa(
          par = rep(0, n_theta),
          fn = current_objective_sequential_dcorICA,
          lower = rep(-pi, n_theta),
          upper = rep(pi, n_theta),
          control = control,
          Z = Z_current,
          transform = transform,
          alpha = alpha,
          block_size = block_size
        )
      } else {
        stop("Invalid optimizer chosen.")
      }
      
      thetas[current_index:(current_index + n_theta - 1L)] <- res$par
      
      if (optimizer == "cobyla"){
        sweep_log[[i]] <- list(
          step = i,
          par = res$par,
          objective = res$value,
          convergence = res$convergence,
          message = if (!is.null(res$message)) res$message else NA_character_
        )
      } else if (optimizer == "bobyqa"){
        sweep_log[[i]] <- list(
          step = i,
          par = res$par,
          objective = res$fval,
          convergence = res$ierr,
          message = if (!is.null(res$msg)) res$msg else NA_character_
        )
      }
      
      Q <- diag(ncol(Z_current))
      for (k in seq_len(n_theta)) {
        Q <- givens_rotation(
          thetas[current_index + k - 1L],
          p = ncol(Z_current),
          index = c(1L + k, 1L)
        ) %*% Q
      }
      
      Z_current <- tcrossprod(Z_current, Q)[, -1L, drop = FALSE]
      current_index <- current_index + n_theta
    }
    
    W_new <- theta_to_W(thetas)[inverse_perm, inverse_perm, drop = FALSE]
    W_total <- W_new %*% W_total
    Z <- tcrossprod(Z, W_new)
    
    convergence_log[[sweep_idx]] <- list(
      sweep = sweep_idx,
      permutation = perm,
      steps = sweep_log
    )
  }
  
  list(
    W_rot = W_total,
    convergence = convergence_log
  )
}

# Check the number of permutation resamples
.check_B <- function(B) {
  if (!is.numeric(B) || length(B) != 1L || is.na(B) || B < 1 || B != as.integer(B)) {
    stop("'B' must be a single positive integer.", call. = FALSE)
  }
  invisible(TRUE)
}

# With K(k/m), m = length(lags), a lag k >= m gets weight zero. This is expected
# for the largest lag of 1:m, but otherwise selected lags are silently ignored.
.warn_kernel_lags <- function(lags) {
  m <- length(lags)
  if (m < 2L || any(lags > m)) {
    warning(
      "Lags larger than the number of selected lags, and a single selected lag, ",
      "get kernel weight zero and do not contribute to the test statistic. ",
      "Consider lags = 1:m with m >= 2 or type = \"BP\".",
      call. = FALSE
    )
  }
  invisible(TRUE)
}

# Check the scale values used to standardize the data
.check_scale_positive <- function(scale) {
  if (anyNA(scale) || any(scale <= 0)) {
    stop(
      "All values of `scale` must be positive. The default MAD is zero if more ",
      "than half of the observations are tied; in this case supply `scale`.",
      call. = FALSE
    )
  }
  invisible(TRUE)
}

# Internal infix helper
`%||%` <- function(x, y) {
  if (is.null(x)) y else x
}

# Coerce vector/matrix input to numeric matrix
.as_matrix_input <- function(X) {
  if (is.vector(X)) {
    X <- matrix(X, ncol = 1L)
  } else {
    X <- as.matrix(X)
  }
  X
}

.biloop_transform_univariate_scaled <- function(x, c1 = 4, c2 = 4) {
  z <- tanh(as.numeric(x) / c1)
  out <- .twocircles_map(z)
  out[, 1] <- c2 * out[, 1]
  out
}

.twocircles_map <- function(x) {
  x <- as.numeric(x)
  xax <- numeric(length(x))
  neg <- x < 0
  
  xax[neg] <- -cos(2 * pi * (x[neg] - 0.5)) - 1
  xax[!neg] <-  cos(2 * pi * (x[!neg] + 0.5)) + 1
  yax <- sin(2 * pi * x)
  
  cbind(xax, yax)
}