diff options
| -rw-r--r-- | DESCRIPTION | 27 | ||||
| -rw-r--r-- | MD5 | 18 | ||||
| -rw-r--r-- | NAMESPACE | 25 | ||||
| -rw-r--r-- | R/check_inputs.R | 64 | ||||
| -rw-r--r-- | R/cv_loss.R | 17 | ||||
| -rw-r--r-- | R/plot_print_methods.R | 26 | ||||
| -rw-r--r-- | R/rmad.R | 31 | ||||
| -rw-r--r-- | R/rsc.R | 42 | ||||
| -rw-r--r-- | R/rsc_cv.R | 183 | ||||
| -rwxr-xr-x | R/zzz.R | 4 | ||||
| -rw-r--r-- | inst/CITATION | 40 | ||||
| -rwxr-xr-x | man/plot.cv_rsc.Rd | 84 | ||||
| -rwxr-xr-x | man/rmad.Rd | 136 | ||||
| -rwxr-xr-x | man/rsc.Rd | 118 | ||||
| -rwxr-xr-x | man/rsc_cv.Rd | 228 | ||||
| -rw-r--r-- | src/Makevars | 9 | ||||
| -rw-r--r-- | src/cormaddp.f90 | 82 | ||||
| -rw-r--r-- | src/init.c | 22 | ||||
| -rw-r--r-- | src/selectiongeneral.f95 | 1081 |
19 files changed, 2237 insertions, 0 deletions
diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..36f47f6 --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,27 @@ +Package: RSC +Type: Package +Title: Robust and Sparse Correlation Matrix +Description: Performs robust and sparse correlation matrix estimation. Robustness is achieved based on a simple robust pairwise correlation estimator, while sparsity is obtained based on thresholding. The optimal thresholding is tuned via cross-validation. See Serra, Coretto, Fratello, and Tagliaferri (2018) <doi:10.1093/bioinformatics/btx642>. +Authors@R: c(person("Luca", "Coraggio", + role = c("cre", "aut"), email = "luca.coraggio@unina.it"), + person("Pietro", "Coretto", + role = c("aut"), email = "pcoretto@unisa.it"), + person("Angela", "Serra", + role = c("aut"), email = "angela.serra@tuni.fi"), + person("Roberto", "Tagliaferri", + role = c("ctb"), email = "robtag@unisa.it")) +Author: Luca Coraggio [cre, aut], + Pietro Coretto [aut], + Angela Serra [aut], + Roberto Tagliaferri [ctb] +Maintainer: Luca Coraggio <luca.coraggio@unina.it> +NeedsCompilation: yes +Imports: stats, graphics, Matrix, methods, parallel, foreach, + doParallel, utils +License: GPL (>= 2) +LazyData: TRUE +Version: 1.0 +Date: 2020-06-29 +Packaged: 2020-07-02 07:30:56 UTC; pietro +Repository: CRAN +Date/Publication: 2020-07-04 10:50:03 UTC @@ -0,0 +1,18 @@ +80c587bf8be31b19549b965aa44f1d50 *DESCRIPTION +1f64d30a870f2eae4855e83010aaaa4e *NAMESPACE +0b87a473118069181446fb007472c042 *R/check_inputs.R +1df1ecba75175b1665069fc88419cf02 *R/cv_loss.R +1ada6bc39f15ea58f34a5d46173141fa *R/plot_print_methods.R +6dcfdaa8b26686b70794ba2c9d6147ca *R/rmad.R +f26fdd51a2378e500e6ed2eaeda8721a *R/rsc.R +556a9a10b85bfad178103f3cc5fdc0ff *R/rsc_cv.R +607039efb30fc06f123290e3dcbce1d4 *R/zzz.R +ee184b1708ad4ef4735070065fe6737c *inst/CITATION +fffcd92fdbda643ea59f12298cc04287 *man/plot.cv_rsc.Rd +57435eb3ed73e1af9435073eeaf57928 *man/rmad.Rd +18e1373418dce90153d0727d19f5fa9e *man/rsc.Rd +2ab80b72b7090b2514c9cf4b3fbe006d *man/rsc_cv.Rd +24cec7b34e7c1931f74d1f25cad05a83 *src/Makevars +f135feaeda1c3e803c74729a94fcc416 *src/cormaddp.f90 +9ded1c66a2fd469e79d7d1d25d0bd604 *src/init.c +62f9244cb578a6720b664233da503e32 *src/selectiongeneral.f95 diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..ffd8478 --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,25 @@ +## import all functions from the following pkgs +import("stats", "graphics", "Matrix") + +## selective imports from pkgs +importFrom("methods", "as") +importFrom("foreach", "foreach", "%dopar%") +importFrom("parallel","detectCores") +importFrom("doParallel", "registerDoParallel", "stopImplicitCluster") +importFrom("utils", "citHeader", "citEntry") + +## export objects in /src +useDynLib(RSC, .registration = TRUE, .fixes = "F_") + +## export objects in /R +export(rmad) +export(rsc_cv) +export(rsc) + + +## export S3 methods +S3method(plot, rsc_cv) +S3method(print, rsc_cv) + +## ## exports every object that doesn't start with a dot +## exportPattern("^[^\\.]") 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) + } +} @@ -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) +} @@ -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() +} diff --git a/inst/CITATION b/inst/CITATION new file mode 100644 index 0000000..8690fa8 --- /dev/null +++ b/inst/CITATION @@ -0,0 +1,40 @@ +citHeader("To cite package 'RSC' in publications use:") + +year <- sub("-.*", "", meta$Date) +note <- sprintf("R package version %s", meta$Version) + +citEntry(entry="Manual", + title = "RSC: robust and sparse correlation matrix estimation", + author = personList(person(given="Luca", family="Coraggio"), + person(given="Pietro", family="Coretto"), + person(given="Angela", family="Serra"), + person(given="Roberto", family="Tagliaferri")), + year = year, + note = note, + textVersion = + paste("Coraggio, L., Coretto, P., Serra, A. and Tagliaferri, R. (", year,"). + RSC: robust and sparse correlation matrix estimation. ", + note, "url: https://CRAN.R-project.org/package=rsc", + sep='') + ) + + + +citEntry(entry="Article", + title = "Robust and sparse correlation matrix estimation for the analysis + of high-dimensional genomics data", + author = personList(person(given="Angela", family="Serra"), + person(given="Pietro", family="Coretto"), + person(given="Michele", family="Fratello"), + person(given="Roberto", family="Tagliaferri") + ), + year = 2018, + journal = "Bioinformatics", + volume = 34, + number = 4, + pages = "625--634", + doi = "10.1093/bioinformatics/btx642", + textVersion = "Serra, A., Coretto, P., Fratello, M., and Tagliaferri, R. (2018). Robust and sparse correlation matrix estimation for the analysis of high-dimensional genomics data. Bioinformatics, 34(4), 625-634." + ) + + diff --git a/man/plot.cv_rsc.Rd b/man/plot.cv_rsc.Rd new file mode 100755 index 0000000..0d7b2c2 --- /dev/null +++ b/man/plot.cv_rsc.Rd @@ -0,0 +1,84 @@ +\name{plot.rsc_cv} +\alias{plot.rsc_cv} + +\title{ + Plot method for rsc_cv objects +} + +\description{ + Plot the cross-validation estimates of the Frobenius loss. +} + + +\usage{ + \method{plot}{rsc_cv}(x, \dots) +} + + + + +\arguments{ + \item{x}{ + Output from \code{\link{rsc_cv}}, that is an S3 object of class \code{"rsc_cv"}. + } + \item{\dots}{ + additional arguments passed to \code{\link[graphics]{plot.default}}. + } +} + + +\value{ + Plot the Frobenius loss estimated via cross-validation (y-axis) vs + threshold values (x-axis). The dotted blue line represents the average + expected normalized Frobenius loss, while the vertical segments + around the average are \emph{1-standard-error} error bars + (see \emph{Details} in \code{\link{rsc_cv}}. + + The vertical dashed red line identifies the minimum of the average + loss, that is the optimal threshold flagged as \code{"minimum"}. The + vertical dashed green line identifies the optimal selection flagged + as \code{"minimum1se"} in the output of \code{\link{rsc_cv}} (see + \emph{Details} in \code{\link{rsc_cv}}). +} + + + + +\section{References}{ + Serra, A., Coretto, P., Fratello, M., and Tagliaferri, R. (2018). + Robust and sparsecorrelation matrix estimation for the analysis of + high-dimensional genomics data. \emph{Bioinformatics}, 34(4), + 625-634. doi:10.1093/bioinformatics/btx642 +} + + +\seealso{ + \code{\link{rsc_cv}} +} + + + +\examples{ +\donttest{ +## simulate a random sample from a multivariate Cauchy distribution +## note: example in high-dimension are obtained increasing p +set.seed(1) +n <- 100 # sample size +p <- 10 # dimension +dat <- matrix(rt(n*p, df = 1), nrow = n, ncol = p) +colnames(dat) <- paste0("Var", 1:p) + + +## perform 10-fold cross-validation repeated R=10 times +## note: for multi-core machines experiment with 'ncores' +set.seed(2) +a <- rsc_cv(x = dat, R = 10, K = 10) +a + +## plot the cross-validation estimates +plot(a) + +## pass additional parameters to graphics::plot +plot(a , cex = 2) +} +} diff --git a/man/rmad.Rd b/man/rmad.Rd new file mode 100755 index 0000000..01247cf --- /dev/null +++ b/man/rmad.Rd @@ -0,0 +1,136 @@ +\name{rmad} + +\alias{rmad} + +\title{RMAD correlation matrix} + +\description{ + Compute the RMAD robust correlation matrix proposed in Serra et + al. (2018) based on the robust correlation coefficient proposed in + Pasman and Shevlyakov (1987). +} + + +\usage{ + rmad(x , y = NULL, na.rm = FALSE , even.correction = FALSE) +} + + +\arguments{ + \item{x}{ + A numeric vector, a matrix or a data.frame. If \code{x} is a matrix + or a data.frame, rows of \code{x} correspond to sample units + and columns correspond to variables. If \code{x} is a numerical + vector, and \code{y} is not \code{NULL}, the RMAD correlation + coefficient between \code{x} and \code{y} is computed. Categorical + variables are not allowed. + } + \item{y}{ + A numerical vector if not \code{NULL}. If both \code{x} and \code{y} + are numerical vectors, the RMAD correlation coefficient between + \code{x} and \code{y} is computed. + } + \item{na.rm}{ + A logical value, if \code{TRUE} sample observation + containing \code{NA} values are excluded (see \emph{Details}). + } + \item{even.correction}{ + A logical value, if \code{TRUE} a correction + for the calculation of the medians is applied to reduce the bias + when the number of samples even (see \emph{Details}). + } +} + + +\details{ + The \code{rmad} function computes the correlation matrix based on the + pairwise robust correlation coefficient of Pasman and Shevlyakov + (1987). This correlation coefficient is based on repeated median + calculations for all pairs of variables. This is a computational + intensive task when the number of variables (that is \code{ncol(x)}) + is large. + + The software is optimized for large dimensional data sets, the median + is approximated as the central observation obtained based on the + \emph{introselect} sorting algorithm of Musser (1997) implemented in + Fortran 95 language. For small samples this may be a crude + approximation, however, it makes the computational cost feasible for + high-dimensional data sets. With the option \code{even.correction + = TRUE} a correction is applied to reduce the bias for data sets with + an even number of samples. Although \code{even.correction = TRUE} + has a small computational cost for each pair of variables, it is + suggested to use the default \code{even.correction = FALSE} for large + dimensional data sets. + + The function can handle a data matrix with missing values (\code{NA} + records). If \code{na.rm = TRUE} then missing values are handled by + casewise deletion (and if there are no complete cases, an error is + returned). In practice, if \code{na.rm = TRUE} all rows of + \code{x} that contain at least an \code{NA} are removed. + + Since the software is optimized to work with high-dimensional data sets, + the output RMAD matrix is packed into a storage efficient format + using the \code{"dspMatrix"} S4 class from the \code{\link{Matrix}} + package. The latter is specifically designed for dense real symmetric + matrices. A sparse correlation matrix can be obtained applying + thresholding using the \code{\link{rsc_cv}} and \code{\link{rsc}}. +} + + + +\value{ + \item{If \code{x} is a matrix or a data.frame}{ + Returns a correlation matrix of class \code{"dspMatrix"} (S4 class object) + as defined in the \code{\link{Matrix}} package. + } + \item{If \code{x} and \code{y} are numerical vectors}{ + Returns a numerical value, that is the RMAD correlation coefficient + between \code{x} and \code{y}. + } +} + + + +\section{References}{ + Musser, D. R. (1997). Introspective sorting and selection algorithms. + \emph{Software: Practice and Experience}, 27(8), 983-993. + + Pasman,V. and Shevlyakov,G. (1987). Robust methods of estimation of + correlation coefficient. \emph{Automation Remote Control}, 48, 332-340. + + Serra, A., Coretto, P., Fratello, M., and Tagliaferri, R. (2018). + Robust and sparsecorrelation matrix estimation for the analysis of + high-dimensional genomics data. \emph{Bioinformatics}, 34(4), 625-634. + doi: 10.1093/bioinformatics/btx642 +} + + + +\seealso{ + \code{rsc_cv}, \code{rsc} +} + + + + + + + +\examples{ +## simulate a random sample from a multivariate Cauchy distribution +set.seed(1) +n <- 100 # sample size +p <- 7 # dimension +dat <- matrix(rt(n*p, df = 1), nrow = n, ncol = p) +colnames(dat) <- paste0("Var", 1:p) + + +## compute the rmad correlation coefficient between dat[,1] and dat[,2] +a <- rmad(x = dat[,1], y = dat[,2]) + + +## compute the RMAD correlaiton matrix +b <- rmad(x = dat) +b +} + diff --git a/man/rsc.Rd b/man/rsc.Rd new file mode 100755 index 0000000..cd252a6 --- /dev/null +++ b/man/rsc.Rd @@ -0,0 +1,118 @@ +\name{rsc} + +\alias{rsc} + +\title{Robust and Sparse Correlation Matrix Estimator} + +\description{ + Compute the Robust and Sparse Correlation Matrix (RSC) estimator + proposed in Serra et al. (2018). +} + + +\usage{ + rsc(cv, threshold = "minimum") +} + + +\arguments{ + \item{cv}{ + An S3 object of class \code{"rsc_cv"} (see \code{\link{rsc_cv}}). + } + \item{threshold}{ + Threshold parameter to compute the RSC estimate. This + is a numeric value taken onto the interval (0,1), or it is + equal to \code{"minimum"} or \code{"minimum1se"} for selecting the + optimal threshold according to the selection performed in + \code{\link{rsc_cv}}. + } +} + + +\details{ + The setting \code{threshold = "minimum"} or \code{threshold = + "minimum1se"} applies thresholding according to the criteria + discussed in the \emph{Details} section in \code{\link{rsc_cv}}. + When \code{cv} is obtained using \code{\link{rsc_cv}} with + \code{cv.type = "random"}, the default settings for \code{\link{rsc}} + implements exactly the RSC estimator proposed in Serra et al., + (2018). + + Although \code{threshold = "minimum"} is the default choice, in + high-dimensional situations \code{threshold = "minimum1se"} usually + provides a more parsimonious representation of the correlation + structure. Since the underlying RMAD matrix is passed through the + \code{cv} input, any other hand-tuned threshold to the RMAD matrix + can be applied without significant additional computational + costs. The latter can be done setting \code{threshold} to any value + onto the (0,1) interval. + + The software is optimized to handle high-dimensional data sets, + therefore, the output RSC matrix is packed into a storage efficient + sparse format using the \code{"dsCMatrix"} S4 class from the + \code{\link{Matrix}} package. The latter is specifically designed for + sparse real symmetric matrices. +} + + + +\value{ + Returns a sparse correlaiton matrix of class \code{"dsCMatrix"} + (S4 class object) as defined in the \code{\link{Matrix}} package. +} + + + + +\section{References}{ + Serra, A., Coretto, P., Fratello, M., and Tagliaferri, R. (2018). + Robust and sparsecorrelation matrix estimation for the analysis of + high-dimensional genomics data. \emph{Bioinformatics}, 34(4), + 625-634. doi:10.1093/bioinformatics/btx642 +} + + + +\seealso{ + \code{\link{rsc_cv}} +} + + + + + + + +\examples{ +\donttest{ +## simulate a random sample from a multivariate Cauchy distribution +## note: example in high-dimension are obtained increasing p +set.seed(1) +n <- 100 # sample size +p <- 10 # dimension +dat <- matrix(rt(n*p, df = 1), nrow = n, ncol = p) +colnames(dat) <- paste0("Var", 1:p) + + +## perform 10-fold cross-validation repeated R=10 times +## note: for multi-core machines experiment with 'ncores' +set.seed(2) +a <- rsc_cv(x = dat, R = 10, K = 10) +a + +## obtain the RSC matrix with "minimum" flagged solution +b <- rsc(cv = a, threshold = "minimum") +b + +## obtain the RSC matrix with "minimum1se" flagged solution +d <- rsc(cv = a, threshold = "minimum1se") +d + +## since the object 'a' stores the RMAD underlying estimator, we can +## apply thresholding at any level without re-estimating the RMAD +## matrix +e <- rsc(cv = a, threshold = 0.5) +e +} +} + diff --git a/man/rsc_cv.Rd b/man/rsc_cv.Rd new file mode 100755 index 0000000..3cd2da6 --- /dev/null +++ b/man/rsc_cv.Rd @@ -0,0 +1,228 @@ +\name{rsc_cv} + +\alias{rsc_cv} + +\title{Optimal threshold selection for the RSC estimator} + +\description{ + Perform cross-validation to select an adaptive optimal threshold for + the RSC estimator proposed in Serra et al. (2018). +} + + +\usage{ + rsc_cv(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) +} + + +\arguments{ + \item{x}{ + A matrix or a data.frame. Rows of \code{x} correspond to sample units + and columns correspond to variables. Categorical variables are not + allowed. + } + \item{cv.type}{ + A character string indicating the cross-validation algorithm. Possible + values are \code{"kfold"} for repeated K-fold cross-validation, and + \code{"random"} for random cross-validation (see \emph{Details}). + } + \item{R}{ + An integer corresponding to the number of repeated foldings when + \code{cv.type = "kfold"}. When \code{cv.type = "random"} \code{R} + defines the number of random splits (see \emph{Details}). + } + \item{K}{ + An integer corresponding to the number of \emph{folds} in K-fold + cross-validation. Therefore this argument is not relevant when + \code{cv.type = "random"}. + } + \item{threshold}{ + A sequence of reals taken onto the interval (0,1) defining the + threshold values at which the loss is estimated. + } + \item{even.correction}{ + A logical value. It sets the parameter \code{even.correction} in + each of the underlying RMAD computations (see \emph{Details} in + \code{\link{rmad}}). + } + \item{na.rm}{ + A logical value, it defines the treatment of missing values in + each of the underlying RMAD computations (see \emph{Details}). + } + \item{ncores}{ + An integer value defining the number of cores used for parallel + computing. When \code{ncores=NULL} (default), the number \code{r} of + available cores is detected, and \code{(r-1)} of them are used + (see \emph{Details}). + } + \item{monitor}{ + A logical value. If \code{TRUE} progress messages are + printed on screen. + } +} + + +\details{ + The \code{rsc_cv} function performs cross-validation to estimate the + expected Frobenius loss proposed in Bickel and Levina (2008). The + original contribution of Bickel and Levina (2008), and its extension + in Serra et al. (2018), is based on a random + cross-validation algorithm where the training/test size depends on + the sample size \emph{n}. The latter is implemented selecting + \code{cv.type = "ramdom"}, and fixing an appropriate number \code{R} of random + train/test splits. \code{R} should be as large as possible, but + in practice this impacts the computing time strongly for + high-dimensional data sets. + + Although Serra et al. (2018) showed that the random cross-validation + of Bickel and Levina (2008) works well for the RSC estimator, + subsequent experiments suggested that repeated K-fold cross-validation + on average produces better results. Repeated K-fold cross-validation + is implemented with the default \code{cv.type = "kfold"}. In this case + \code{K} defines the number of \emph{folds}, while \code{R} defines + the number of times that the K-fold cross-validation is repeated with + \code{R} independent shuffles of the original data. Selecting + \code{R=1} and \code{K=10} one performs the standard 10-fold + cross-validation. Ten replicates (\code{R=10}) of the K-fold + cross-validation are generally sufficient to obtain reasonable + estimates of the underlying loss, but for extremely high-dimensional + data \code{R} may be varied to speed up calculations. + + On multi-core hardware the cross-validation is executed in parallel + setting \code{ncores}. The parallelism is implemented on the + total number of data splits, that is \code{R} for the random + cross-validation, and \code{R*K} for the repeated K-fold + cross-validation. The software is optimized so that generally the + total computing time scales almost linearly with the number of + available computer cores (\code{ncores}). + + For both the random and the K-fold cross-validation it is computed the + normalized version of the expected squared Frobenius loss proposed in + Bickel and Levina (2008). The normalization is such + that the squared Frobenius norm of the identity matrix equals to 1 + whatever is its dimension. + + Two optimal threshold selection types are reported with flags (see + \emph{Value} section below): \code{"minimum"} and + \code{"minimum1se"}. The flag \code{"minimum"} denotes the threshold + value that minimizes the average loss. The flag \code{"minimum1se"} + implements the so called + \emph{1-SE rule}: this is the maximum threshold value such that the + corresponding average loss is within \emph{1-standard-error} with + respect to the threshold that minimizes the average loss + (that is the one corresponding to the \code{"minimum"} flag). + + Since unbiased standard errors for the K-fold cross-validation are + impossible to compute (see Bengio and Grandvalet, 2004), when + \code{cv.type="kfold"} the reported standard errors have to be + considered as a downward biased approximation. +} + + + +\value{ + An S3 object of class \code{'cv_rsc'} with the following components: + \item{rmadvec}{ + A vector containing the lower triangle of the underlying RMAD + matrix. + } + \item{varnames}{ + A character vector if variable names are available for the input + data set \code{x}. Otherwise this is \code{NULL}. + } + \item{loss}{ + A data.frame reporting cross-validation estimates. Columns of + \code{loss} are as follows: \code{loss$Threshold} is the threshold value; + \code{loss$Average} is averaged loss; \code{loss$SE} is the standard error + for the average loss; \code{loss$Flag="minimum"} denotes the threshold + achieving the minimum average loss; \code{loss$Flag = "*"} denotes threshold + values such that the average loss is within \emph{1-standard-error} + with respect to the \code{"minimum"} solution. + } + \item{minimum}{ + A numeric value. This is the minimum of the average loss. This + corresponds to the flag \code{"minimum"} in the loss component + above (see \emph{Details}). + } + \item{minimum1se}{ + A numeric value. This is the largest threshold such that the + corresponding \code{flag = "*"}. In practice this selects the + optimal threshold based on the \emph{1-SE rule} discussed in the + \emph{Details} Section above. + } +} + + + + +\section{References}{ + Bengio, Y., and Grandvalet, Y. (2004). No unbiased estimator of the + variance of k-fold cross-validation. \emph{Journal of Machine Learning + Research}, 5(Sep), 1089-1105. + + Bickel, P. J., and Levina, E. (2008). Covariance regularization by + thresholding. The \emph{Annals of Statistics}, 36(6), 2577-2604. + doi:10.1214/08-AOS600 + + Serra, A., Coretto, P., Fratello, M., and Tagliaferri, R. (2018). + Robust and sparsecorrelation matrix estimation for the analysis of + high-dimensional genomics data. \emph{Bioinformatics}, 34(4), + 625-634. doi:10.1093/bioinformatics/btx642 +} + + + +\seealso{ + \code{rsc}, \code{plot.rsc_cv} +} + + + + + + + +\examples{ +\donttest{ +## simulate a random sample from a multivariate Cauchy distribution +## note: example in high-dimension are obtained increasing p +set.seed(1) +n <- 100 # sample size +p <- 10 # dimension +dat <- matrix(rt(n*p, df = 1), nrow = n, ncol = p) +colnames(dat) <- paste0("Var", 1:p) + + +## perform 10-fold cross-validation repeated R=10 times +## note: for multi-core machines experiment with 'ncores' +set.seed(2) +a <- rsc_cv(x = dat, R = 10, K = 10) +a + + +## threshold selection: note that here, knowing the sampling designs, +## we would like to threshold any correlation larger than zero in +## absolute value +## +a$minimum ## "minimum" flagged solution +a$minimum1se ## "minimum1se" flagged solution + +## plot the cross-validation estimates +plot(a) + +## to obtain the RSC matrix we pass 'a' to the rsc() function +b <- rsc(cv = a, threshold = "minimum") +b + +d <- rsc(cv = a, threshold = "minimum1se") +d + +## since the object 'a' stores the RMAD underlying estimator, we can +## apply thresholding at any level without re-estimating the RMAD +## matrix +e <- rsc(cv = a, threshold = 0.5) +e +} +} + diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..8d466ea --- /dev/null +++ b/src/Makevars @@ -0,0 +1,9 @@ +all: $(SHLIB) + $(MAKE) $(SHLIB) + rm -f *.o *.mod + +cormaddp.o: selectiongeneral.o + $(FC) -c -fPIC cormaddp.f90 -o cormaddp.o + +selectiongeneral.o: selectiongeneral.f95 + $(FC) -c -fPIC selectiongeneral.f95 -o selectiongeneral.o diff --git a/src/cormaddp.f90 b/src/cormaddp.f90 new file mode 100644 index 0000000..f9a0f69 --- /dev/null +++ b/src/cormaddp.f90 @@ -0,0 +1,82 @@ + subroutine cormadvecdp(matrix,nrow,ncol,res,ressize,correcteven) + + use selectionalgo + + implicit none + + !note kind dp is defined in selectionalgo + integer, intent(in) :: nrow, ncol,ressize, correcteven !.Fortran R function will not deal with deferred size arrays + real(kind=dp), dimension(nrow,ncol), intent(inout) :: matrix + real(kind=dp), dimension(ressize), intent(out) :: res + + integer :: i, j, n, p + real(kind=dp) :: med,mad, fresh + real(kind=dp), dimension(:), allocatable :: U, V, A, B + + + p=ubound(matrix,2) + n=ubound(matrix,1) + allocate(U(n),V(n)) + + if (correcteven==1) then + do i=1,p-3,3 + U=matrix(:,i) + med=iselect(U,evencorrection=.true.) + matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.)) + U=matrix(:,i+1) + med=iselect(U,evencorrection=.true.) + matrix(:,i+1)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.)) + U=matrix(:,i+2) + med=iselect(U,evencorrection=.true.) + matrix(:,i+2)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.)) + end do + do i=i,p + U=matrix(:,i) + med=iselect(U,evencorrection=.true.) + matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med),evencorrection=.true.)) + end do + !unrolled loops + n=1 + do i=1,p-1 + do j=i+1,p + U=matrix(:,i)+matrix(:,j) + V=-matrix(:,i)+matrix(:,j) + mad=(cost*iselect(abs(U-iselect(U,evencorrection=.true.)),evencorrection=.true.))**2 + med=(cost*iselect(abs(V-iselect(V,evencorrection=.true.)),evencorrection=.true.))**2 + res(n)=(mad-med)/(mad+med) + n=n+1 + end do + end do + else + do i=1,p-3,3 + U=matrix(:,i) + med=iselect(U) + matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med))) + U=matrix(:,i+1) + med=iselect(U) + matrix(:,i+1)=(U-med)/(sqrt2*cost*iselect(abs(U-med))) + U=matrix(:,i+2) + med=iselect(U) + matrix(:,i+2)=(U-med)/(sqrt2*cost*iselect(abs(U-med))) + end do + do i=i,p + U=matrix(:,i) + med=iselect(U) + matrix(:,i)=(U-med)/(sqrt2*cost*iselect(abs(U-med))) + end do + !unrolled loops + n=1 + do i=1,p-1 + do j=i+1,p + U=matrix(:,i)+matrix(:,j) + V=-matrix(:,i)+matrix(:,j) + mad=(cost*iselect(abs(U-iselect(U))))**2 + med=(cost*iselect(abs(V-iselect(V))))**2 + res(n)=(mad-med)/(mad+med) + n=n+1 + end do + end do + end if + + + end subroutine cormadvecdp diff --git a/src/init.c b/src/init.c new file mode 100644 index 0000000..66209f7 --- /dev/null +++ b/src/init.c @@ -0,0 +1,22 @@ +#include <stdlib.h> +//#include <stdio.h> +#include <R.h> +#include <Rinternals.h> +#include <Rmath.h> +#include <Rdefines.h> +#include <R_ext/RS.h> +#include <R_ext/Rdynload.h> + +extern void F77_NAME(cormadvecdp)(void *, void *, void *, void *, void *, void *); + +static const R_FortranMethodDef FortMethods[] = { + {"cormadvecdp", (DL_FUNC) &F77_NAME(cormadvecdp), 6}, + {NULL, NULL, 0} +}; + +void +R_init_RSC(DllInfo *dll) +{ + R_registerRoutines(dll, NULL, NULL, FortMethods, NULL); + R_useDynamicSymbols(dll, FALSE); +} diff --git a/src/selectiongeneral.f95 b/src/selectiongeneral.f95 new file mode 100644 index 0000000..1078261 --- /dev/null +++ b/src/selectiongeneral.f95 @@ -0,0 +1,1081 @@ +module selectionalgo + use, intrinsic :: iso_fortran_env + + implicit none + + real, parameter :: cost=1.4826, sqrt2=sqrt(2.) + integer, parameter :: sp = real32, dp=real64 + + interface qselect + procedure quickselectsp, quickselectscalarsp,quickselectdp, quickselectscalardp + end interface qselect + + interface iselect + procedure introselectsp, introselectscalarsp, introselectdp, introselectscalardp + end interface iselect + + contains + + !!!!!!! SINGLE PRECISION + !from here there are the selectionalgo and selectionalgoimpl condensed in one file + recursive subroutine quickselectrecursivesp(invector,tovector,k,output) + real(kind=sp), dimension(:), allocatable, target, intent(inout) :: invector + real(kind=sp), intent(inout), pointer, contiguous :: tovector(:) + integer, intent(inout) :: k + real(kind=sp), intent(out) :: output + real(kind=sp) :: swapper, pvt + integer :: i, r, subst + real(kind=sp), pointer :: a,b,c + + r=size(tovector) + + select case(r) + case(1) + output=tovector(1) + return + case(2) + select case(k) + case(1) + output=minval(tovector) + case(2) + output=maxval(tovector) + end select + return + case(3) + select case(k) + case(1) + output=minval(tovector) + case(3) + output=maxval(tovector) + case(2) + i=r/2+1 + output=sum(tovector)-minval(tovector)-maxval(tovector) + end select + return + case default + i=r/2+1 + + !pivoting section (pivot of 3) + a=>tovector(1) + b=>tovector(i) + c=>tovector(r) + pvt=a+b+c + swapper=max(a,b,c) + b=swapper + pvt=pvt-swapper + swapper=min(a,b,c) + a=swapper + pvt=pvt-swapper + c=pvt + + !one pass section + subst=1 + b=>tovector(subst) + do i=1,(r-1-3),3 + a => tovector(i) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + a => tovector(i+1) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + a => tovector(i+2) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + end do + do i=i,(r-1) + a => tovector(i) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + end do + tovector(r)=tovector(subst) + tovector(subst)=pvt + + !decide next step + select case(subst-k) + case(0) + output=tovector(k) + case(1:) + tovector => tovector(1:subst-1) + call quickselectrecursivesp(invector,tovector,k, output) + case(:-1) + k=k-subst + tovector => tovector(subst+1:) + call quickselectrecursivesp(invector,tovector,k, output) + end select + end select + end subroutine quickselectrecursivesp + + recursive subroutine introselectrecursivesp(invector,tovector,k,output) + real(kind=sp), dimension(:), allocatable, target, intent(inout) :: invector + real(kind=sp), intent(inout), pointer, contiguous :: tovector(:) + integer, intent(inout) :: k + real(kind=sp), intent(out) :: output + real(kind=sp) :: swapper, pvt + integer, save :: switch + integer :: i, r, subst, r_min, switch_after + real(kind=sp), pointer :: a,b,c + + r=ubound(tovector,1) + r_min = 3000 ! at least this bigger to switch to median of medians + switch_after = 5 ! at least fails this amount of time to reorder + ! one third of the vector to consider med of med + + select case(r) + case(1) + output=tovector(1) + switch=0 + return + case(2) + select case(k) + case(1) + output=minval(tovector) + case(2) + output=maxval(tovector) + end select + switch=0 + return + case(3) + select case(k) + case(1) + output=minval(tovector) + case(3) + output=maxval(tovector) + case(2) + i=r/2+1 + output=sum(tovector)-minval(tovector)-maxval(tovector) + end select + switch=0 + return + case default + i=r/2+1 + + !pivoting section (pivot of 3) + a=>tovector(1) + b=>tovector(i) + c=>tovector(r) + pvt=a+b+c + swapper=max(a,b,c) + b=swapper + pvt=pvt-swapper + swapper=min(a,b,c) + a=swapper + pvt=pvt-swapper + c=pvt + + !one pass section + subst=1 + b=>tovector(subst) + do i=1,(r-1-3),3 + a => tovector(i) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + a => tovector(i+1) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + a => tovector(i+2) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + end do + do i=i,(r-1) + a => tovector(i) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + end do + tovector(r)=tovector(subst) + tovector(subst)=pvt + + !decide next step + select case(subst-k) + case(0) + output=tovector(k) + switch=0 + case(1:) + tovector => tovector(1:subst-1) + if ((subst-1)>r*2/3) then + switch=switch+1 + if (switch==switch_after .and. r>r_min) then + !reset switch and call medofmedsp + switch=0 + call medofmedsp(invector,tovector,k, output) + else + call introselectrecursivesp(invector,tovector,k, output) + end if + else + switch=0 + call introselectrecursivesp(invector,tovector,k, output) + end if + case(:-1) + k=k-subst + tovector => tovector(subst+1:) + if ((r-subst)>r*2/3) then + switch=switch+1 + if (switch==switch_after .and. r>r_min) then + !reset switch and call medofmedsp + switch=0 + call medofmedsp(invector,tovector,k, output) + else + call introselectrecursivesp(invector,tovector,k, output) + end if + else + switch=0 + call introselectrecursivesp(invector,tovector,k, output) + end if + end select + end select + end subroutine introselectrecursivesp + + recursive subroutine medofmedsp(invector,tovector,k,output) + real(kind=sp), dimension(:), allocatable, target, intent(inout) :: invector + real(kind=sp), intent(inout), pointer, contiguous :: tovector(:) + real(kind=sp), intent(out) :: output + integer, intent(inout) :: k + + real(kind=sp), pointer, contiguous :: subsec(:) + real(kind=sp), dimension(:),allocatable :: meds + integer :: i,subst,r + real(kind=sp) :: pvt, swapper + real(kind=sp), pointer :: a,b + + r=ubound(tovector,1) + select case(mod(r,5)) + case(0) + allocate(meds(r/5)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(i/5+1)=pvt + end do + case(1) + allocate(meds(r/5+1)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(i/5+1)=pvt + end do + meds(r/5+1)=tovector(i) + case(2) + allocate(meds(r/5+1)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(i/5+1)=pvt + end do + subsec=>tovector(i:) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + meds(r/5+1)=subsec(2) + case(3) + allocate(meds(r/5+1)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(i/5+1)=pvt + end do + subsec=>tovector(i:) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + meds(r/5+1)=pvt + case(4) + allocate(meds(r/5+1)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(i/5+1)=pvt + end do + subsec=>tovector(i:) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(r/5+1)=pvt + end select + + !set pivot + pvt=qselect(meds) + do i=1,r-3,3 + if (tovector(i)==pvt) then + tovector(i)=tovector(r) + tovector(r)=pvt + exit + else if (tovector(i+1)==pvt) then + tovector(i+1)=tovector(r) + tovector(r)=pvt + exit + else if (tovector(i+2)==pvt) then + tovector(i+2)=tovector(r) + tovector(r)=pvt + exit + else + cycle + end if + end do + do i=i,r + if (tovector(i)==pvt) then + tovector(i)=tovector(r) + tovector(r)=pvt + exit + else + cycle + end if + end do + + call introselectrecursivesp(invector,tovector,k, output) + end subroutine medofmedsp + + function introselectsp(invector,ord,evencorrection) + real(kind=sp), intent(in), dimension(:) :: invector + logical, optional, intent(in) :: evencorrection + integer, intent(in), optional :: ord + real(kind=sp) :: introselectsp + + real(kind=sp), dimension(:), allocatable, target :: vector + real(kind=sp), pointer, contiguous :: tovector(:) => null() + integer :: k1, k + + k1=size(invector) + allocate(vector(k1)) + vector=invector + tovector => vector + + k1=ubound(invector,1) + if(present(ord)) then + k=ord + else + k=k1/2+1 + end if + + if(present(evencorrection)) then + if (k/=k1/2+1 .or. mod(k1,2)/=0) then + !write(*,'(a)') 'Warning: evencorrection argument ignored.' + k1=0 + else + k1=0 + if (evencorrection) k1=k-1 + end if + else + k1=0 + end if + + call introselectrecursivesp(vector,tovector,k, introselectsp) + + if (k1/=0) then + block + real(kind=sp) :: temp + temp=introselectsp + tovector=>vector + call introselectrecursivesp(vector,tovector,k1, introselectsp) + introselectsp=(introselectsp+temp)/2.0 + end block + end if + end function introselectsp + + function quickselectsp(invector,ord,evencorrection) + real(kind=sp), intent(in), dimension(:) :: invector + logical, optional, intent(in) :: evencorrection + integer, intent(in), optional :: ord + real(kind=sp) :: quickselectsp + + real(kind=sp), dimension(:), allocatable, target :: vector + real(kind=sp), pointer, contiguous :: tovector(:) => null() + integer :: k1, k + + k1=size(invector) + allocate(vector(k1)) + vector=invector + tovector => vector + + k1=ubound(invector,1) + if(present(ord)) then + k=ord + else + k=k1/2+1 + end if + + if(present(evencorrection)) then + if (k/=k1/2+1 .or. mod(k1,2)/=0) then + !write(*,'(a)') 'Warning: evencorrection argument ignored.' + k1=0 + else + k1=0 + if (evencorrection) k1=k-1 + end if + else + k1=0 + end if + + call quickselectrecursivesp(vector,tovector,k, quickselectsp) + + if (k1/=0) then + block + real(kind=sp) :: temp + temp=quickselectsp + tovector=>vector(1:k1+1) + call quickselectrecursivesp(vector,tovector,k1, quickselectsp) + quickselectsp=(quickselectsp+temp)/2.0 + end block + end if + end function quickselectsp + + !deal with scalar input + function introselectscalarsp(invector,ord, evencorrection) + real(kind=sp), intent(in) :: invector + integer, intent(in), optional :: ord + logical, optional, intent(in) :: evencorrection + real(kind=sp) :: introselectscalarsp + + introselectscalarsp=invector + end function introselectscalarsp + + function quickselectscalarsp(invector,ord, evencorrection) + real(kind=sp), intent(in) :: invector + integer, intent(in), optional :: ord + logical, optional, intent(in) :: evencorrection + real(kind=sp) :: quickselectscalarsp + + quickselectscalarsp=invector + end function quickselectscalarsp + + !!!!!!! DOUBLE PRECISION + !from here there are the selectionalgo and selectionalgoimpl condensed in one file + recursive subroutine quickselectrecursivedp(invector,tovector,k,output) + real(kind=dp), dimension(:), allocatable, target, intent(inout) :: invector + real(kind=dp), intent(inout), pointer, contiguous :: tovector(:) + integer, intent(inout) :: k + real(kind=dp), intent(out) :: output + real(kind=dp) :: swapper, pvt + integer :: i, r, subst + real(kind=dp), pointer :: a,b,c + + r=size(tovector) + + select case(r) + case(1) + output=tovector(1) + return + case(2) + select case(k) + case(1) + output=minval(tovector) + case(2) + output=maxval(tovector) + end select + return + case(3) + select case(k) + case(1) + output=minval(tovector) + case(3) + output=maxval(tovector) + case(2) + i=r/2+1 + output=sum(tovector)-minval(tovector)-maxval(tovector) + end select + return + case default + i=r/2+1 + + !pivoting section (pivot of 3) + a=>tovector(1) + b=>tovector(i) + c=>tovector(r) + pvt=a+b+c + swapper=max(a,b,c) + b=swapper + pvt=pvt-swapper + swapper=min(a,b,c) + a=swapper + pvt=pvt-swapper + c=pvt + + !one pass section + subst=1 + b=>tovector(subst) + do i=1,(r-1-3),3 + a => tovector(i) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + a => tovector(i+1) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + a => tovector(i+2) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + end do + do i=i,(r-1) + a => tovector(i) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + end do + tovector(r)=tovector(subst) + tovector(subst)=pvt + + !decide next step + select case(subst-k) + case(0) + output=tovector(k) + case(1:) + tovector => tovector(1:subst-1) + call quickselectrecursivedp(invector,tovector,k, output) + case(:-1) + k=k-subst + tovector => tovector(subst+1:) + call quickselectrecursivedp(invector,tovector,k, output) + end select + end select + end subroutine quickselectrecursivedp + + recursive subroutine introselectrecursivedp(invector,tovector,k,output) + real(kind=dp), dimension(:), allocatable, target, intent(inout) :: invector + real(kind=dp), intent(inout), pointer, contiguous :: tovector(:) + integer, intent(inout) :: k + real(kind=dp), intent(out) :: output + real(kind=dp) :: swapper, pvt + integer, save :: switch + integer :: i, r, subst, r_min, switch_after + real(kind=dp), pointer :: a,b,c + + r=ubound(tovector,1) + r_min = 3000 ! at least this bigger to switch to median of medians + switch_after = 5 ! at least fails this amount of time to reorder + ! one third of the vector to consider med of med + + select case(r) + case(1) + output=tovector(1) + switch=0 + return + case(2) + select case(k) + case(1) + output=minval(tovector) + case(2) + output=maxval(tovector) + end select + switch=0 + return + case(3) + select case(k) + case(1) + output=minval(tovector) + case(3) + output=maxval(tovector) + case(2) + i=r/2+1 + output=sum(tovector)-minval(tovector)-maxval(tovector) + end select + switch=0 + return + case default + i=r/2+1 + + !pivoting section (pivot of 3) + a=>tovector(1) + b=>tovector(i) + c=>tovector(r) + pvt=a+b+c + swapper=max(a,b,c) + b=swapper + pvt=pvt-swapper + swapper=min(a,b,c) + a=swapper + pvt=pvt-swapper + c=pvt + + !one pass section + subst=1 + b=>tovector(subst) + do i=1,(r-1-3),3 + a => tovector(i) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + a => tovector(i+1) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + a => tovector(i+2) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + end do + do i=i,(r-1) + a => tovector(i) + if (a<pvt) then + swapper=b + b=a + a=swapper + subst=subst+1 + b=>tovector(subst) + end if + end do + tovector(r)=tovector(subst) + tovector(subst)=pvt + + !decide next step + select case(subst-k) + case(0) + output=tovector(k) + switch=0 + case(1:) + tovector => tovector(1:subst-1) + if ((subst-1)>r*2/3) then + switch=switch+1 + if (switch==switch_after .and. r>r_min) then + !reset switch and call medofmeddp + switch=0 + call medofmeddp(invector,tovector,k, output) + else + call introselectrecursivedp(invector,tovector,k, output) + end if + else + switch=0 + call introselectrecursivedp(invector,tovector,k, output) + end if + case(:-1) + k=k-subst + tovector => tovector(subst+1:) + if ((r-subst)>r*2/3) then + switch=switch+1 + if (switch==switch_after .and. r>r_min) then + !reset switch and call medofmeddp + switch=0 + call medofmeddp(invector,tovector,k, output) + else + call introselectrecursivedp(invector,tovector,k, output) + end if + else + switch=0 + call introselectrecursivedp(invector,tovector,k, output) + end if + end select + end select + end subroutine introselectrecursivedp + + recursive subroutine medofmeddp(invector,tovector,k,output) + real(kind=dp), dimension(:), allocatable, target, intent(inout) :: invector + real(kind=dp), intent(inout), pointer, contiguous :: tovector(:) + real(kind=dp), intent(out) :: output + integer, intent(inout) :: k + + real(kind=dp), pointer, contiguous :: subsec(:) + real(kind=dp), dimension(:),allocatable :: meds + integer :: i,subst,r + real(kind=dp) :: pvt, swapper + real(kind=dp), pointer :: a,b + + r=ubound(tovector,1) + select case(mod(r,5)) + case(0) + allocate(meds(r/5)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(i/5+1)=pvt + end do + case(1) + allocate(meds(r/5+1)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(i/5+1)=pvt + end do + meds(r/5+1)=tovector(i) + case(2) + allocate(meds(r/5+1)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(i/5+1)=pvt + end do + subsec=>tovector(i:) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + meds(r/5+1)=subsec(2) + case(3) + allocate(meds(r/5+1)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(i/5+1)=pvt + end do + subsec=>tovector(i:) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + meds(r/5+1)=pvt + case(4) + allocate(meds(r/5+1)) + do i=1,r-4,5 + subsec=>tovector(i:i+4) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(i/5+1)=pvt + end do + subsec=>tovector(i:) + subst=minloc(subsec,1) + pvt=subsec(subst) + subsec(subst)=subsec(1) + subsec(1)=pvt + subst=minloc(subsec(2:),1)+1 + pvt=subsec(subst) + subsec(subst)=subsec(2) + subsec(2)=pvt + subst=minloc(subsec(3:),1)+2 + pvt=subsec(subst) + subsec(subst)=subsec(3) + subsec(3)=pvt + meds(r/5+1)=pvt + end select + + !set pivot + pvt=qselect(meds) + do i=1,r-3,3 + if (tovector(i)==pvt) then + tovector(i)=tovector(r) + tovector(r)=pvt + exit + else if (tovector(i+1)==pvt) then + tovector(i+1)=tovector(r) + tovector(r)=pvt + exit + else if (tovector(i+2)==pvt) then + tovector(i+2)=tovector(r) + tovector(r)=pvt + exit + else + cycle + end if + end do + do i=i,r + if (tovector(i)==pvt) then + tovector(i)=tovector(r) + tovector(r)=pvt + exit + else + cycle + end if + end do + + call introselectrecursivedp(invector,tovector,k, output) + end subroutine medofmeddp + + function introselectdp(invector,ord,evencorrection) + real(kind=dp), intent(in), dimension(:) :: invector + logical, optional, intent(in) :: evencorrection + integer, intent(in), optional :: ord + real(kind=dp) :: introselectdp + + real(kind=dp), dimension(:), allocatable, target :: vector + real(kind=dp), pointer, contiguous :: tovector(:) => null() + integer :: k1, k + + k1=size(invector) + allocate(vector(k1)) + vector=invector + tovector => vector + + k1=ubound(invector,1) + if(present(ord)) then + k=ord + else + k=k1/2+1 + end if + + if(present(evencorrection)) then + if (k/=k1/2+1 .or. mod(k1,2)/=0) then + !write(*,'(a)') 'Warning: evencorrection argument ignored.' + k1=0 + else + k1=0 + if (evencorrection) k1=k-1 + end if + else + k1=0 + end if + + call introselectrecursivedp(vector,tovector,k, introselectdp) + + if (k1/=0) then + block + real(kind=dp) :: temp + temp=introselectdp + !tovector=>vector(1:k1+1) + tovector=>vector + call introselectrecursivedp(vector,tovector,k1, introselectdp) + introselectdp=(introselectdp+temp)/2.0 + end block + end if + end function introselectdp + + function quickselectdp(invector,ord,evencorrection) + real(kind=dp), intent(in), dimension(:) :: invector + logical, optional, intent(in) :: evencorrection + integer, intent(in), optional :: ord + real(kind=dp) :: quickselectdp + + real(kind=dp), dimension(:), allocatable, target :: vector + real(kind=dp), pointer, contiguous :: tovector(:) => null() + integer :: k1, k + + k1=size(invector) + allocate(vector(k1)) + vector=invector + tovector => vector + + k1=ubound(invector,1) + if(present(ord)) then + k=ord + else + k=k1/2+1 + end if + + if(present(evencorrection)) then + if (k/=k1/2+1 .or. mod(k1,2)/=0) then + !write(*,'(a)') 'Warning: evencorrection argument ignored.' + k1=0 + else + k1=0 + if (evencorrection) k1=k-1 + end if + else + k1=0 + end if + + call quickselectrecursivedp(vector,tovector,k, quickselectdp) + + if (k1/=0) then + block + real(kind=dp) :: temp + temp=quickselectdp + tovector=>vector(1:k1+1) + call quickselectrecursivedp(vector,tovector,k1, quickselectdp) + quickselectdp=(quickselectdp+temp)/2.0 + end block + end if + end function quickselectdp + + !deal with scalar input + function introselectscalardp(invector,ord, evencorrection) + real(kind=dp), intent(in) :: invector + integer, intent(in), optional :: ord + logical, optional, intent(in) :: evencorrection + real(kind=dp) :: introselectscalardp + + introselectscalardp=invector + end function introselectscalardp + + function quickselectscalardp(invector,ord, evencorrection) + real(kind=dp), intent(in) :: invector + integer, intent(in), optional :: ord + logical, optional, intent(in) :: evencorrection + real(kind=dp) :: quickselectscalardp + + quickselectscalardp=invector + end function quickselectscalardp + + +end module selectionalgo |