version 1.0

7 files changed, 367 insertions, 0 deletions

diff --git a/R/check_inputs.R b/R/check_inputs.R
new file mode 100644
index 0000000..4cbcf0a
--- /dev/null
+++ b/R/check_inputs.R
@@ -0,0 +1,64 @@
+.check_input_data_matrix <- function(x, y, na.rm) {
+  if (is.null(y)) {
+    if (!is.array(x) & !is.data.frame(x) & !is.matrix(x)) {
+      stop("\"x\" must be a numeric matrix or any other array-type object that can be converted or a \"matrix\" object with with ncol>=2.\n\n")
+    }
+    if (is.vector(x)) {
+      stop("\"x\" must be a numeric matrix or any other array-type object that can be converted or a \"matrix\" object.  with with ncol>=2.\n\n")
+    }
+    if (!is.matrix(x)) {
+      x <- data.matrix(x)
+    }
+    if (!is.numeric(x)) {
+      stop("\"x\" must be numeric.")
+    }
+    if (nrow(x) < 2 | ncol(x) < 2) {
+      stop("nrow(x)>=2 and ncol(xa)>=2 are required\n\n")
+    }
+    is_na_data <- is.na(x)
+    if (any(is_na_data)) {
+      if (na.rm == FALSE) {
+        stop("\"x\" contains NA records. You may want to filter NAs by setting \"na.rm=TRUE\" (see documentation for more details).\n\n")
+      }
+      else {
+        idx_na <- which(rowSums(is_na_data) >= 1)
+        x <- x[-idx_na, , drop = FALSE]
+        if (nrow(x) < 2) {
+          stop("nrow(x)<2 after NA removal.\n\n")
+        }
+      }
+    }
+    if (any(!is.finite(x))) {
+      stop("\"x\" contains Inf values\n\n")
+    }
+  }
+  else {
+    if (!is.vector(x) | !is.vector(y)) {
+      stop("If \"y\" is given, \"x\" and  \"y\" must be both numeric.\n\n")
+    }
+    if (!is.numeric(x) | !is.numeric(y)) {
+      stop("\"x\" and \"y\" must be numeric.\n\n")
+    }
+    if (length(x) != length(y)) {
+      stop("\"x\" and \"y\" have different length.\n\n")
+    }
+    x <- cbind(x, y, deparse.level = 0)
+    is_na_data <- is.na(x)
+    if (any(is_na_data)) {
+      if (na.rm == FALSE) {
+        stop("\"x\" or \"y\" contains NA records. You may want to filter NAs by setting \"na.rm=TRUE\" (see documentation for more details).\n\n")
+      }
+      else {
+        idx_na <- which(rowSums(is_na_data) >= 1)
+        x <- x[-idx_na, , drop = FALSE]
+        if (nrow(x) < 2) {
+          stop("length(x)<2 and/or length(y)<2 after NA removal.\n\n")
+        }
+      }
+    }
+    if (any(!is.finite(x))) {
+      stop("\"x\" and/or \"y\"  contains Inf values\n\n")
+    }
+  }
+  return(x)
+}
diff --git a/R/cv_loss.R b/R/cv_loss.R
new file mode 100644
index 0000000..ea46007
--- /dev/null
+++ b/R/cv_loss.R
@@ -0,0 +1,17 @@
+.cv_loss <- function(idx, dat, evencorrection, threshold, grid.length, p, nc) {
+  res <- numeric(nc)
+  n1 <- as.integer(sum(idx))
+  C1 <- .Fortran("cormadvecdp", matrix = dat[idx, ], nrow = n1, ncol = p, res = res, 
+    ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res
+  n2 <- as.integer(sum(!idx))
+  C2 <- .Fortran("cormadvecdp", matrix = dat[!idx, ], nrow = n2, ncol = p, res = res, 
+    ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res
+  ans <- rep(0, times = grid.length)
+  for (h in 1:grid.length) {
+    C1[abs(C1) < threshold[h]] <- 0
+    ans[h] <- sum(2 * {
+      C1 - C2
+    }^2)/p
+  }
+  return(ans)
+}
diff --git a/R/plot_print_methods.R b/R/plot_print_methods.R
new file mode 100644
index 0000000..ca163cb
--- /dev/null
+++ b/R/plot_print_methods.R
@@ -0,0 +1,26 @@
+print.rsc_cv <- function(x, ...) {
+  cat("\n")
+  cat("================================================\n")
+  cat("   Expected Normalized Squared Frobenius Loss   \n")
+  cat("================================================\n")
+  print(x$loss)
+  cat("================================================\n")
+  cat("\n")
+}
+plot.rsc_cv <- function(x, ...) {
+  tstar <- which(x$loss$Flag == "minimum")
+  hstar <- x$loss$Threshold[tstar]
+  inf_loss <- x$loss$Average - x$loss$SE
+  sup_loss <- x$loss$Average + x$loss$SE
+  a <- inf_loss[tstar]
+  b <- sup_loss[tstar]
+  hstar1se <- max(x$loss$Threshold[which(x$loss$Flag == "*")])
+  Ylim <- range(c(inf_loss, sup_loss))
+  plot(x$loss$Threshold, x$loss$Average, t = "b", ylim = Ylim, pch = 20, col = "#0052A5", 
+    lwd = 2, main = "RSC Optimal Threshold Selection", xlab = "Threshold", ylab = "Average loss", 
+    ...)
+  arrows(x$loss$Threshold, inf_loss, x$loss$Threshold, sup_loss, length = 0.05, 
+    angle = 90, code = 3, col = "#0052A5")
+  abline(v = hstar1se, col = "#31A853", lty = 2, lwd = 2)
+  abline(v = hstar, col = "#E0162B", lty = 2, lwd = 2)
+}
diff --git a/R/rmad.R b/R/rmad.R
new file mode 100644
index 0000000..b8e9d41
--- /dev/null
+++ b/R/rmad.R
@@ -0,0 +1,31 @@
+rmad <- function(x, y = NULL, na.rm = FALSE, even.correction = FALSE) {
+  dat <- .check_input_data_matrix(x = x, y = y, na.rm = na.rm)
+  colnames_original <- colnames(dat)
+  storage.mode(dat) <- "double"
+  n <- as.integer(nrow(dat))
+  p <- as.integer(ncol(dat))
+  nc <- as.integer({
+    p^2 - p
+  }/2)
+  if (even.correction) {
+    evencorrection <- 1L
+  }
+  else {
+    evencorrection <- 0L
+  }
+  u <- .Fortran("cormadvecdp", matrix = dat, nrow = n, ncol = p, res = numeric(nc), 
+    ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res
+  if (!is.null(y)) {
+    return(u)
+  }
+  else {
+    R <- Matrix(1, nrow = p, ncol = p, sparse = FALSE)
+    R[lower.tri(R, diag = FALSE)] <- u
+    R <- forceSymmetric(R, uplo = "L")
+    R <- as(R, "dspMatrix")
+    if (!is.null(colnames_original)) {
+      dimnames(R)[[1]] <- dimnames(R)[[2]] <- colnames_original
+    }
+    return(R)
+  }
+}
diff --git a/R/rsc.R b/R/rsc.R
new file mode 100644
index 0000000..6cf4e15
--- /dev/null
+++ b/R/rsc.R
@@ -0,0 +1,42 @@
+rsc <- function(cv, threshold = "minimum") {
+  if (class(cv) == "rsc_cv") {
+    if (is.numeric(threshold)) {
+      if (length(threshold) > 1) {
+        stop("if a specific value for 'threshold' is chosen, this must be a single numeric value in (0,1)")
+      }
+      else if (threshold <= 0 | threshold >= 1) {
+        stop("if a specific value for 'threshold' is chosen, this must be a single numeric value in (0,1)")
+      }
+    }
+    else {
+      if ({
+        threshold != "minimum"
+      } & {
+        threshold != "minimum1se"
+      }) {
+        stop("'threshold' must be one of the following: 'minimum', 'minimum1se', a numeric value in (0,1).")
+      }
+      if (threshold == "minimum") {
+        threshold <- cv$minimum
+      }
+      else if (threshold == "minimum1se") {
+        threshold <- cv$minimum1se
+      }
+    }
+    cv$rmadvec[abs(cv$rmadvec) < threshold] <- 0
+    nc <- length(cv$rmadvec)
+    p <- {
+      1 + sqrt(1 + 8 * nc)
+    }/2
+    R <- Matrix(1, nrow = p, ncol = p, sparse = TRUE)
+    R[lower.tri(R, diag = FALSE)] <- cv$rmadvec
+    R <- forceSymmetric(R, uplo = "L")
+    if (!is.null(cv$varnames)) {
+      dimnames(R)[[1]] <- dimnames(R)[[2]] <- cv$varnames
+    }
+  }
+  else {
+    stop("'cv' must be a an object of class 'rsc_cv' obtained from 'rsc::rsc_cv'")
+  }
+  return(R)
+}
diff --git a/R/rsc_cv.R b/R/rsc_cv.R
new file mode 100644
index 0000000..4153322
--- /dev/null
+++ b/R/rsc_cv.R
@@ -0,0 +1,183 @@
+rsc_cv <- function(x, cv.type = "kfold", R = 10, K = 10, threshold = seq(0.05, 0.95, 
+  by = 0.025), even.correction = FALSE, na.rm = FALSE, ncores = NULL, monitor = TRUE) {
+  dat <- .check_input_data_matrix(x = x, y = NULL, na.rm = na.rm)
+  colnames_original <- colnames(dat)
+  storage.mode(dat) <- "double"
+  n <- as.integer(nrow(dat))
+  p <- as.integer(ncol(dat))
+  nc <- as.integer({
+    p^2 - p
+  }/2)
+  if ({
+    cv.type != "random"
+  } & {
+    cv.type != "kfold"
+  }) {
+    stop("\"cv.type\" must be either \"random\" (default)  or \"kfold\"")
+  }
+  if (!is.numeric(R)) {
+    stop("\"R\" must be an integer > 1")
+  }
+  else if (R < 1) {
+    stop("\"R\" must be an integer > 1")
+  }
+  if (!is.numeric(K)) {
+    stop("\"K\" must be an integer > 1")
+  }
+  else if (R < 1) {
+    stop("\"K\" must be an integer > 1")
+  }
+  if (length(threshold) == 1) {
+    if (threshold <= 0 | threshold >= 1) {
+      stop("\"threshold\" value does not belong to the interval (0,1).")
+    }
+  }
+  else if (length(threshold) > 1) {
+    if (any(threshold < 0) | any(threshold >= 1)) {
+      stop("Some of the \"threshold\" values do not belong to the interval (0,1).")
+    }
+    grid.length <- length(threshold)
+  }
+  if (even.correction) {
+    evencorrection <- 1L
+  }
+  else {
+    evencorrection <- 0L
+  }
+  if (is.null(ncores)) {
+    DetectedCores <- detectCores()
+    if (DetectedCores <= 2) {
+      ncores <- 1
+    }
+    else {
+      ncores <- {
+        DetectedCores - 1
+      }
+    }
+  }
+  else {
+    ncores <- as.integer(ncores)
+    if (ncores <= 0) {
+      stop("\"ncores\" must be an integer larger or equal to 1.")
+    }
+  }
+  if (monitor) {
+    cat("\n")
+    message("Computing the RMAD matrix")
+    t0 <- Sys.time()
+  }
+  rmad_vec <- .Fortran("cormadvecdp", matrix = dat, nrow = n, ncol = p, res = numeric(nc), 
+    ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res
+  if (monitor) {
+    t1 <- Sys.time()
+    dt01 <- difftime(t1, t0, units = "auto")
+    message("* RMAD computing time:...... ", round(dt01, 2), " [", attributes(dt01)$units, 
+      "]")
+  }
+  if (cv.type == "random") {
+    if (monitor) {
+      cat("\n")
+      message("Performing cross-validation")
+      t_hat <- round(1.2 * {
+        {
+          dt01 * R * 2
+        }/ncores
+      }, 2)
+      message("* predicted end time (worst case):...... ", Sys.time() + t_hat)
+    }
+    n1 <- n - floor(n/log(n))
+    IDX <- array(FALSE, dim = c(R, n))
+    for (r in 1:R) {
+      IDX[r, ][sample(1:n, size = n1, replace = FALSE)] <- TRUE
+    }
+    registerDoParallel(ncores)
+    U <- foreach(r = 1:R) %dopar% {
+      .cv_loss(idx = IDX[r, ], dat = dat, evencorrection = evencorrection, 
+        threshold = threshold, grid.length = grid.length, p = p, nc = nc)
+    }
+    stopImplicitCluster()
+    LOSS <- array(0, dim = c(R, grid.length))
+    for (r in 1:R) {
+      LOSS[r, ] <- U[[r]]
+    }
+    avg_loss <- apply(LOSS, 2, mean)
+    se_loss <- apply(LOSS, 2, sd)/sqrt(R)
+  }
+  if (cv.type == "kfold") {
+    if (monitor) {
+      cat("\n")
+      message("Performing cross-validation")
+      t_hat <- round(1.2 * {
+        {
+          dt01 * R * K * 2
+        }/ncores
+      }, 2)
+      message("* predicted end time (worst case):...... ", Sys.time() + t_hat)
+    }
+    idx_fold <- cut(1:n, breaks = K, labels = FALSE)
+    IDX <- array(TRUE, dim = c(R * K, n))
+    row_count <- 1L
+    for (r in 1:R) {
+      idx_fold_shuffle <- sample(idx_fold, size = n, replace = FALSE)
+      for (k in 1:K) {
+        IDX[row_count, ][idx_fold_shuffle == k] <- FALSE
+        row_count <- 1L + row_count
+      }
+    }
+    registerDoParallel(ncores)
+    U <- foreach(r = 1:{
+      R * K
+    }) %dopar% {
+      .cv_loss(idx = IDX[r, ], dat = dat, evencorrection = evencorrection, 
+        threshold = threshold, grid.length = grid.length, p = p, nc = nc)
+    }
+    stopImplicitCluster()
+    if (R == 1) {
+      LOSS <- array(0, dim = c(K, grid.length))
+      for (k in 1:K) {
+        LOSS[k, ] <- U[[k]]
+      }
+    }
+    else {
+      LOSS <- array(0, dim = c(K, grid.length, R))
+      dimnames(LOSS)[[3]] <- paste0("r", 1:R)
+      dimnames(LOSS)[[1]] <- paste0("k", 1:K)
+      row_count <- 1L
+      for (r in 1:R) {
+        for (k in 1:K) {
+          LOSS[k, , r] <- U[[row_count]]
+          row_count <- 1L + row_count
+        }
+      }
+    }
+    avg_loss_r <- sd_loss_r <- matrix(0, nrow = R, ncol = grid.length)
+    for (r in 1:R) {
+      avg_loss_r[r, ] <- apply(LOSS[, , r], 2, mean)
+      sd_loss_r[r, ] <- apply(LOSS[, , r], 2, sd)/sqrt(K)
+    }
+    avg_loss <- apply(avg_loss_r, 2, mean)
+    se_loss <- apply(sd_loss_r, 2, mean)
+    if (monitor) {
+      t2 <- Sys.time()
+      dt02 <- difftime(t2, t0, units = "auto")
+    }
+  }
+  if (monitor) {
+    message("* finished on:.......................... ", Sys.time())
+  }
+  tstar <- which.min(avg_loss)
+  flags <- rep("", grid.length)
+  a <- avg_loss[tstar] - se_loss[tstar]
+  b <- avg_loss[tstar] + se_loss[tstar]
+  flags[{
+    avg_loss >= a
+  } & {
+    avg_loss <= b
+  }] <- "*"
+  flags[tstar] <- "minimum"
+  res <- data.frame(Threshold = threshold, Average = avg_loss, SE = se_loss, Flag = flags)
+  ans <- list(rmadvec = rmad_vec, varnames = colnames_original, loss = res, minimum = threshold[tstar], 
+    minimum1se = max(threshold[avg_loss >= a & avg_loss <= b]))
+  class(ans) <- "rsc_cv"
+  return(ans)
+}
diff --git a/R/zzz.R b/R/zzz.R
new file mode 100755
index 0000000..3c33499
--- /dev/null
+++ b/R/zzz.R
@@ -0,0 +1,4 @@
+.onAttach <- function(lib, pkg) {
+  packageStartupMessage("\nRSC: robust and sparse correlation matrix estimation\n     Type 'citation(\"RSC\")' for citing this package\n")
+  invisible()
+}