From ca17d57dce048f57e03241f6120d539ec70d785a Mon Sep 17 00:00:00 2001 From: Luca Coraggio Date: Sun, 17 Oct 2021 19:00:08 +0000 Subject: version 2.0 --- DESCRIPTION | 9 +- MD5 | 29 ++- NAMESPACE | 2 +- NEWS | 10 + R/check_inputs.R | 139 +++++----- R/cv_loss.R | 28 +- R/plot_print_methods.R | 63 +++-- R/rmad.R | 42 ++- R/rsc.R | 89 ++++--- R/rsc_cv.R | 226 ++++++++++++---- R/zzz.R | 9 +- man/rmad.Rd | 11 +- src/Makevars | 2 + src/RSCdefines.h | 11 + src/cormad_ptr.c | 180 +++++++++++++ src/cormaddp.f90 | 642 ---------------------------------------------- src/init.c | 74 ++++-- src/selection_algos_ptr.c | 154 +++++++++++ 18 files changed, 849 insertions(+), 871 deletions(-) create mode 100644 src/Makevars create mode 100644 src/RSCdefines.h create mode 100644 src/cormad_ptr.c delete mode 100644 src/cormaddp.f90 create mode 100644 src/selection_algos_ptr.c diff --git a/DESCRIPTION b/DESCRIPTION index 3f611ee..dcafcfb 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -19,9 +19,8 @@ NeedsCompilation: yes Imports: stats, graphics, Matrix, methods, parallel, foreach, doParallel, utils License: GPL (>= 2) -LazyData: TRUE -Version: 1.2 -Date: 2020-07-24 -Packaged: 2020-07-24 07:47:23 UTC; pietro +Version: 2.0 +Date: 2021-10-14 +Packaged: 2021-10-14 14:44:39 UTC; luco Repository: CRAN -Date/Publication: 2020-07-24 14:00:07 UTC +Date/Publication: 2021-10-17 20:00:08 UTC diff --git a/MD5 b/MD5 index 7aaaa51..c43967d 100644 --- a/MD5 +++ b/MD5 @@ -1,17 +1,20 @@ -7a9c09c1b58001a6ef17502cccb6b9ee *DESCRIPTION -1f64d30a870f2eae4855e83010aaaa4e *NAMESPACE -1eae13396706a60fc95f5ff4b751765b *NEWS -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 +f8d0f0d74876f46c0818b142ffa180aa *DESCRIPTION +9e6c9e63967c10a2dad469679534e3b2 *NAMESPACE +12aa4cb77fcbd1096f7c868a332debac *NEWS +852c9f1717e4dbd26ed8c549edc3c588 *R/check_inputs.R +228fb1b443a4c9eaf4fddcdf9c8a31a9 *R/cv_loss.R +8c8c70764c8a9ab72f23ba031e0f330d *R/plot_print_methods.R +720ba075c0b402a232b5528774fa1230 *R/rmad.R +3e98806cb071a0b84028b7e5bb022cc5 *R/rsc.R +98992e4ed337da761f641cb2960f7579 *R/rsc_cv.R +cafa1636fd3af0dd260cbfcb27df210e *R/zzz.R 02a06053793e54ce413a378a05ccc3ce *inst/CITATION 1e6af804927098cfa6c235278682578c *man/plot.cv_rsc.Rd -57435eb3ed73e1af9435073eeaf57928 *man/rmad.Rd +a762f06952b9b9f89cab9826674be146 *man/rmad.Rd 8cc5c9373c7aa8835e9818d0de274ece *man/rsc.Rd 63d3e58e9bdc0da9969feb6ce69d005f *man/rsc_cv.Rd -f45a14eed865277afe99ebcb1bed66f8 *src/cormaddp.f90 -9ded1c66a2fd469e79d7d1d25d0bd604 *src/init.c +95e3011e37d9dde0d75f3a3819b2acd3 *src/Makevars +4aad47542a9a6ce7b9f153dccd1690aa *src/RSCdefines.h +7ef7aff199844c0d42d07684e58e63d1 *src/cormad_ptr.c +b518af5dd0a97f37ab37e4d915dc5303 *src/init.c +cdb3d966f6db89bf54d3fac4b9314bc6 *src/selection_algos_ptr.c diff --git a/NAMESPACE b/NAMESPACE index ffd8478..3d5ba56 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -9,7 +9,7 @@ importFrom("doParallel", "registerDoParallel", "stopImplicitCluster") importFrom("utils", "citHeader", "citEntry") ## export objects in /src -useDynLib(RSC, .registration = TRUE, .fixes = "F_") +useDynLib(RSC, .registration = TRUE, .fixes = "C_") ## export objects in /R export(rmad) diff --git a/NEWS b/NEWS index d7efdcf..f79bd3a 100755 --- a/NEWS +++ b/NEWS @@ -1,3 +1,13 @@ +## 14/10/2021 at 11:51:10 (CEST) +* New version: v2.0. + + Implemented the low-level routines in C. + + Moved from introselect to quickselect (median calculation); + faster on benchmarks with intended input. + + Added support for parallel execution. + +## 07/10/2021 at 11:51:10 (CEST) +* v1.3 does not compile on Solaris (syntax issues) + ## 24/07/2020 at 09:44:39 (CEST) * v1.2 uploaded * fixed compiler issue on Solaris diff --git a/R/check_inputs.R b/R/check_inputs.R index 4cbcf0a..ccc54d4 100644 --- a/R/check_inputs.R +++ b/R/check_inputs.R @@ -1,64 +1,81 @@ -.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") +## Check input data and return a valid matrix object +## +.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 (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(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 (any(!is.finite(x))) { - stop("\"x\" and/or \"y\" contains Inf values\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) -} + + + return(x) + +} ### END function + diff --git a/R/cv_loss.R b/R/cv_loss.R index ea46007..d855362 100644 --- a/R/cv_loss.R +++ b/R/cv_loss.R @@ -1,17 +1,29 @@ -.cv_loss <- function(idx, dat, evencorrection, threshold, grid.length, p, nc) { - res <- numeric(nc) +## Normalized squared Frobenius loss for all threshold values at a given train/test +## set, where +## +## * idx = TRUE for sample points into the train set +## * train set = dat[ idx , ] +## * test set = dat[ !idx , ] +## +.cv_loss <- function(idx, dat, evencorrection, threshold, grid.length, p) { + + + ## compute RMAD on train set n1 <- as.integer(sum(idx)) - C1 <- .Fortran("cormadvecdp", matrix = dat[idx, ], nrow = n1, ncol = p, res = res, - ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res + C1 <- .Call(C_cormad_C, dat[idx, ], n1, p, evencorrection, num.threads = 1) + + ## compute RMAD on test set n2 <- as.integer(sum(!idx)) - C2 <- .Fortran("cormadvecdp", matrix = dat[!idx, ], nrow = n2, ncol = p, res = res, - ressize = nc, correcteven = evencorrection, PACKAGE = "RSC")$res + C2 <- .Call(C_cormad_C, dat[!idx, ], n2, p, correcteven = evencorrection, num.threads = 1) + + ## apply thresholds ans <- rep(0, times = grid.length) for (h in 1:grid.length) { - C1[abs(C1) < threshold[h]] <- 0 + C1[abs(C1) < threshold[h]] <- 0 ## fit on train set ans[h] <- sum(2 * { C1 - C2 - }^2)/p + }^2) / p } + return(ans) } diff --git a/R/plot_print_methods.R b/R/plot_print_methods.R index ca163cb..746179a 100644 --- a/R/plot_print_methods.R +++ b/R/plot_print_methods.R @@ -1,26 +1,41 @@ -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") + +## print method for the crsc_cv class +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) + + + +## ## ## plot method for the crsc_cv class +plot.rsc_cv <- function(x, ...){ + + ## add check object + + 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=.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 index b8e9d41..46e9b35 100644 --- a/R/rmad.R +++ b/R/rmad.R @@ -1,31 +1,47 @@ -rmad <- function(x, y = NULL, na.rm = FALSE, even.correction = FALSE) { +rmad <- function(x, y = NULL, na.rm = FALSE, even.correction = FALSE, num.threads = "half-max") { + + ## check input data 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) + + + ## set even correction if (even.correction) { evencorrection <- 1L - } - else { + } 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) + + ## set number of threads + if (num.threads == "half-max") { + num.threads <- 0L + } else { + storage.mode(num.threads) <- "integer" } - else { + + + ## Call C code + u <- .Call(C_cormad_C, dat, n, p, evencorrection, num.threads) + + + if (!is.null(y)) { ## 2-dimensional + return(u) + } else { ## p-dimensional + + ## assemble the matrix using the lower triangle 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") + + ## attach dimnames if needed if (!is.null(colnames_original)) { dimnames(R)[[1]] <- dimnames(R)[[2]] <- colnames_original } + return(R) - } -} + } ## END if(!is.null(y)){ ## 2-dimensional +} ## END function diff --git a/R/rsc.R b/R/rsc.R index 6cf4e15..8182d1c 100644 --- a/R/rsc.R +++ b/R/rsc.R @@ -1,42 +1,49 @@ -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) +rsc <- function(cv, threshold = "minimum"){ + + ## inputs + ## cv = u ## a class cv_rsc or any other correlation matrix + ## threshold = "minimum" ## "minimum", "minimum1se" or numeric in (0,1) + + if(class(cv) == "rsc_cv"){ + + ## check threshold + 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 + } + } + + + ## threshold the rmadvec + 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") + + ## attach dimnames if needed + 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 index 4153322..660695e 100644 --- a/R/rsc_cv.R +++ b/R/rsc_cv.R @@ -1,79 +1,140 @@ -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) { +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) { + + ## ## ***************************************************************************** + ## ## RSC inputs + ## ## ***************************************************************************** + ## x = X + ## cv.type = "random" ## "kfold" "random" + ## R = 10 ## replicate (for K-fold) or splits for random + ## K = 10 ## folds in kfcv + ## threshold = seq(0.025, 0.975, by = 0.025) + ## opt = "min" , ## "min" "min1se" + ## even.correction = FALSE + ## na.rm = TRUE + ## ncores = 6 + ## monitor = TRUE + ## ## ***************************************************************************** + + + ## check input data 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) + + ## check cv.type if ({ cv.type != "random" } & { cv.type != "kfold" }) { - stop("\"cv.type\" must be either \"random\" (default) or \"kfold\"") + stop('"cv.type" must be either "random" (default) or "kfold"') } + + + ## check R if (!is.numeric(R)) { - stop("\"R\" must be an integer > 1") - } - else if (R < 1) { - stop("\"R\" must be an integer > 1") + stop('"R" must be an integer > 1') + } else if (R < 1) { + stop('"R" must be an integer > 1') } + + + ## check K if (!is.numeric(K)) { - stop("\"K\" must be an integer > 1") - } - else if (R < 1) { - stop("\"K\" must be an integer > 1") + stop('"K" must be an integer > 1') + } else if (R < 1) { + stop('"K" must be an integer > 1') } + + + ## check threshold if (length(threshold) == 1) { if (threshold <= 0 | threshold >= 1) { - stop("\"threshold\" value does not belong to the interval (0,1).") + stop('"threshold" value does not belong to the interval (0,1).') } - } - else if (length(threshold) > 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).") + stop('Some of the "threshold" values do not belong to the interval (0,1).') } grid.length <- length(threshold) } + + + + + + ## set even correction if (even.correction) { evencorrection <- 1L - } - else { + } else { evencorrection <- 0L } + + + + ## set and check ncores if (is.null(ncores)) { DetectedCores <- detectCores() if (DetectedCores <= 2) { ncores <- 1 - } - else { + } else { ncores <- { DetectedCores - 1 } } - } - else { + } else { ncores <- as.integer(ncores) if (ncores <= 0) { - stop("\"ncores\" must be an integer larger or equal to 1.") + stop('"ncores" must be an integer larger or equal to 1.') } } + + + + + + + + + + + ## compute the RMAD vec form 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 + + rmad_vec <- .Call(C_cormad_C, dat, n, p, evencorrection, num.threads = 1) + if (monitor) { t1 <- Sys.time() dt01 <- difftime(t1, t0, units = "auto") - message("* RMAD computing time:...... ", round(dt01, 2), " [", attributes(dt01)$units, - "]") + message( + "* RMAD computing time:...... ", round(dt01, 2), + " [", attributes(dt01)$units, "]" + ) } + + + + + + + + + if (cv.type == "random") { if (monitor) { cat("\n") @@ -81,28 +142,55 @@ rsc_cv <- function(x, cv.type = "kfold", R = 10, K = 10, threshold = seq(0.05, 0 t_hat <- round(1.2 * { { dt01 * R * 2 - }/ncores + } / ncores }, 2) message("* predicted end time (worst case):...... ", Sys.time() + t_hat) } - n1 <- n - floor(n/log(n)) + + + ## IDX[k,i] = TRUE means dat[i, ] is in train set at the k-th split + n1 <- n - floor(n / log(n)) ## train set IDX <- array(FALSE, dim = c(R, n)) for (r in 1:R) { IDX[r, ][sample(1:n, size = n1, replace = FALSE)] <- TRUE } + + ## register parallel backend registerDoParallel(ncores) + + ## parallel computation of losses over splits 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) + .cv_loss( + idx = IDX[r, ], dat = dat, evencorrection = evencorrection, + threshold = threshold, grid.length = grid.length, + p = p + ) } + + ## stop parallel backend stopImplicitCluster() + + ## LOSS[split , threshold] + ## array with normalized squared Frobenius loss LOSS <- array(0, dim = c(R, grid.length)) for (r in 1:R) { LOSS[r, ] <- U[[r]] } + + ## Estimate average loss with std errors avg_loss <- apply(LOSS, 2, mean) - se_loss <- apply(LOSS, 2, sd)/sqrt(R) + se_loss <- apply(LOSS, 2, sd) / sqrt(R) } + + + + + + + + + + if (cv.type == "kfold") { if (monitor) { cat("\n") @@ -110,35 +198,59 @@ rsc_cv <- function(x, cv.type = "kfold", R = 10, K = 10, threshold = seq(0.05, 0 t_hat <- round(1.2 * { { dt01 * R * K * 2 - }/ncores + } / ncores }, 2) message("* predicted end time (worst case):...... ", Sys.time() + t_hat) } + + ## create K deterministic folds idx_fold <- cut(1:n, breaks = K, labels = FALSE) + + ## IDX :: indexes of all train sets + ## IDX[k , i] = TRUE means dat[i, ] is in the train set at some k-th split + ## rows {{r-1}*K + 1 }:{r*K} of IDX correspond to the r-th replica IDX <- array(TRUE, dim = c(R * K, n)) + + ## set initial IDX row counter row_count <- 1L for (r in 1:R) { + ## at each replicate r shuffle the original fold indexes idx_fold_shuffle <- sample(idx_fold, size = n, replace = FALSE) + + ## for each fold shuffle make up the K-fold for (k in 1:K) { IDX[row_count, ][idx_fold_shuffle == k] <- FALSE row_count <- 1L + row_count } } + + + ## register parallel backend registerDoParallel(ncores) + + ## parallel computation of losses over splits 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) + .cv_loss( + idx = IDX[r, ], dat = dat, evencorrection = evencorrection, + threshold = threshold, grid.length = grid.length, + p = p + ) } + + ## stop parallel backend stopImplicitCluster() + + + ## LOSS[fold , threshold, replica] if (R == 1) { LOSS <- array(0, dim = c(K, grid.length)) + ## names ?? for (k in 1:K) { LOSS[k, ] <- U[[k]] } - } - else { + } else { LOSS <- array(0, dim = c(K, grid.length, R)) dimnames(LOSS)[[3]] <- paste0("r", 1:R) dimnames(LOSS)[[1]] <- paste0("k", 1:K) @@ -150,21 +262,33 @@ rsc_cv <- function(x, cv.type = "kfold", R = 10, K = 10, threshold = seq(0.05, 0 } } } + + + ## compute average losses and standard errors 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) + 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") } - } + } ## END if(cv.type == "kfold"){ + + + if (monitor) { message("* finished on:.......................... ", Sys.time()) } + + + ## Organize cross-validation results tstar <- which.min(avg_loss) flags <- rep("", grid.length) a <- avg_loss[tstar] - se_loss[tstar] @@ -175,9 +299,23 @@ rsc_cv <- function(x, cv.type = "kfold", R = 10, K = 10, threshold = seq(0.05, 0 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])) + res <- data.frame( + Threshold = threshold, + Average = avg_loss, + SE = se_loss, + Flag = flags + ) + + ## output list + 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 index 3c33499..ca8c78e 100755 --- a/R/zzz.R +++ b/R/zzz.R @@ -1,4 +1,7 @@ -.onAttach <- function(lib, pkg) { - packageStartupMessage("\nRSC: robust and sparse correlation matrix estimation\n Type 'citation(\"RSC\")' for citing this package\n") - invisible() +.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/man/rmad.Rd b/man/rmad.Rd index 01247cf..26417a0 100755 --- a/man/rmad.Rd +++ b/man/rmad.Rd @@ -12,7 +12,7 @@ \usage{ - rmad(x , y = NULL, na.rm = FALSE , even.correction = FALSE) + rmad(x , y = NULL, na.rm = FALSE , even.correction = FALSE, num.threads = "half-max") } @@ -39,6 +39,9 @@ for the calculation of the medians is applied to reduce the bias when the number of samples even (see \emph{Details}). } + \item{num.threads}{ + An integer value or the string \code{"half-max"} (default), specifying the number of threads for parallel execution (see \emph{Details}). + } } @@ -74,6 +77,12 @@ 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}}. + + \code{rmad} function supports parallel execution. + This is provided via \emph{openmp} (http://www.openmp.org), which must be available; otherwise, falls back to single-core execution. + If \code{num.threads > 0}, function is executed using \code{min(num.threads, max.threads)} threads, where \code{max.threads} is the maximum number of available threads. That is, if positive use the specified number of threads (up to the maximum available). + If \code{num.threads < 0}, function is executed using \code{max(max.threads - num.threads, 1)} threads, i.e. when negative \code{num.threads} indicates the number of threads not to use. + If \code{num.threads == 0} or \code{num.threads == "half-max"}, function is executed using half of the available threads (\code{max(max.threads/2, 1)}). This is the default. } diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..0ee3615 --- /dev/null +++ b/src/Makevars @@ -0,0 +1,2 @@ +PKG_CFLAGS = $(SHLIB_OPENMP_CFLAGS) +PKG_LIBS = $(SHLIB_OPENMP_CFLAGS) diff --git a/src/RSCdefines.h b/src/RSCdefines.h new file mode 100644 index 0000000..75203ac --- /dev/null +++ b/src/RSCdefines.h @@ -0,0 +1,11 @@ +/* Defines quickselect and cormad (double precision) */ + +double quickselect_recursive(double *vector_extract_k, int vec_size, int extract_this_element); + +void cormad(double *matrix, int n_row, int n_col, double *output, int evencorrect); + +#ifdef _OPENMP +void cormad_parallel(double *matrix, int n_row, int n_col, double *output, + int evencorrect, int num_threads); +#endif + diff --git a/src/cormad_ptr.c b/src/cormad_ptr.c new file mode 100644 index 0000000..670a71a --- /dev/null +++ b/src/cormad_ptr.c @@ -0,0 +1,180 @@ +#include "RSCdefines.h" +#include +#include +#ifdef _OPENMP +#include +#endif + +#define SQRT2 (sqrt(2)) +#define CONST 1.4826 + +void cormad(double *matrix, int n_row, int n_col, double *output, + int evencorrect) { + + int k = n_row / 2; // position of the median + int output_size = (n_col - 1) * n_col / 2; + double med, mad; // store variables for medians + double U[n_row], V[n_row]; // help vectors + + /* Transform matrix columns (CORMAD part 1) + * each column has n entries; + * matrix is assumed to be streamed in onevector + */ + int i = 0; + for (int l = 0; l < n_col * n_row; l++) { + if (i == (n_row - 1)) { + U[i] = *matrix; + med = quickselect_recursive(U, n_row, k); + if (evencorrect == 1) { // handle even correction + med = (med + quickselect_recursive(U, n_row, k - 1)) / 2; + } + for (int j = 0; j < n_row; j++) { + U[j] = *(matrix - (n_row - 1) + j) - med; + V[j] = fabs(U[j]); + } + med = quickselect_recursive(V, n_row, k); + if (evencorrect == 1) { // handle even correction + med = (med + quickselect_recursive(V, n_row, k - 1)) / 2; + } + for (int j = 0; j < n_row; j++) // reassign + *(matrix - (n_row - 1) + j) = U[j] / (SQRT2 * CONST * med); + // prepare next iter + matrix++; + i = 0; + } else { + U[i] = *matrix; + // prepare next iter + matrix++; + i++; + } + } + + matrix = matrix - (n_col * n_row); /* reset pointer */ + + /* Operate on columns pairs (CORMAD part 2) */ + double *matrix_2 = matrix; // used to point at second column + int first_col = 0; /* Running first column */ + int second_col = 0; /* Running second column */ + for (int l = 0; l < output_size; l++) { + if (second_col == n_col - 1) { + first_col++; + second_col = first_col; + /* set pointers to columns */ + matrix += n_row; + matrix_2 = matrix; + } + second_col++; + matrix_2 += n_row; + + for (int i = 0; i < n_row; i++) { // auxiliary vectors from matrix + U[i] = *(matrix + i) + *(matrix_2 + i); + V[i] = -*(matrix + i) + *(matrix_2 + i); + } + mad = quickselect_recursive(U, n_row, k); + med = quickselect_recursive(V, n_row, k); + if (evencorrect == 1) { + mad = (mad + quickselect_recursive(U, n_row, k - 1)) / 2; + med = (med + quickselect_recursive(V, n_row, k - 1)) / 2; + } + for (int i = 0; i < n_row; i++) { // reassign for new medians + U[i] = fabs(U[i] - mad); + V[i] = fabs(V[i] - med); + } + mad = quickselect_recursive(U, n_row, k); + med = quickselect_recursive(V, n_row, k); + if (evencorrect == 1) { + mad = (mad + quickselect_recursive(U, n_row, k - 1)) / 2; + med = (med + quickselect_recursive(V, n_row, k - 1)) / 2; + } + mad = pow(CONST * mad, 2); + med = pow(CONST * med, 2); + + // Assign output + *output = (mad - med) / (mad + med); + output++; + } +} + +#ifdef _OPENMP +void cormad_parallel(double *matrix, int n_row, int n_col, double *output, + int evencorrect, int num_threads) { + int k = n_row / 2; // position of the median + int output_size = (n_col - 1) * n_col / 2; + double med, mad; // store variables for medians + double U[n_row], V[n_row]; // help vectors + + /* Transform matrix columns (CORMAD part 1) */ + double *help_matrix = matrix; /* help pointer for matrix */ +#pragma omp parallel for private(med, mad, U, V, help_matrix) \ + num_threads(num_threads) + for (int j = 0; j < n_col; j++) { // iterate on cols + help_matrix = matrix + n_row * j; /* set pointer at beg of col */ + for (int i = 0; i < n_row; i++) + U[i] = *(help_matrix + i); + med = quickselect_recursive(U, n_row, k); + if (evencorrect == 1) { // handle even correction + med = (med + quickselect_recursive(U, n_row, k - 1)) / 2; + } + for (int i = 0; i < n_row; i++) { + U[i] = *(help_matrix + i) - med; + V[i] = fabs(U[i]); + } + med = quickselect_recursive(V, n_row, k); + if (evencorrect == 1) { // handle even correction + med = (med + quickselect_recursive(V, n_row, k - 1)) / 2; + } + for (int i = 0; i < n_row; i++) // reassign + *(help_matrix + i) = U[i] / (SQRT2 * CONST * med); + } + + // int l = 0; // used to iterate over output + /* Operate on columns pairs (CORMAD part 2) */ + help_matrix = matrix; /* help pointers for matrix */ + double *help_matrix_2 = matrix; +#pragma omp parallel for num_threads(num_threads) private( \ + med, mad, U, V, help_matrix, help_matrix_2) + for (int l = 0; l < output_size; l++) { + /* Detrmine columns pairs */ + int col1 = 0, col2 = 0; + int copy_l = l + 1; + for (int last_elem = n_col - 1; last_elem > 0; last_elem--) { + copy_l -= last_elem; + if (copy_l <= 0) { + col2 = (n_col - 1) + copy_l; + break; + } else { + col1++; + } + } + /* Set pointers to columns */ + help_matrix = matrix + n_row * col1; // col1 + help_matrix_2 = matrix + n_row * col2; // col2 + + for (int i = 0; i < n_row; i++) { // auxiliary vectors from matrix + U[i] = *(help_matrix + i) + *(help_matrix_2 + i); + V[i] = -*(help_matrix + i) + *(help_matrix_2 + i); + } + mad = quickselect_recursive(U, n_row, k); + med = quickselect_recursive(V, n_row, k); + if (evencorrect == 1) { + mad = (mad + quickselect_recursive(U, n_row, k - 1)) / 2; + med = (med + quickselect_recursive(V, n_row, k - 1)) / 2; + } + for (int i = 0; i < n_row; i++) { // reassign for new medians + U[i] = fabs(U[i] - mad); + V[i] = fabs(V[i] - med); + } + mad = quickselect_recursive(U, n_row, k); + med = quickselect_recursive(V, n_row, k); + if (evencorrect == 1) { + mad = (mad + quickselect_recursive(U, n_row, k - 1)) / 2; + med = (med + quickselect_recursive(V, n_row, k - 1)) / 2; + } + mad = pow(CONST * mad, 2); + med = pow(CONST * med, 2); + + // Assign output + *(output + l) = (mad - med) / (mad + med); + } +} +#endif diff --git a/src/cormaddp.f90 b/src/cormaddp.f90 deleted file mode 100644 index b0b6516..0000000 --- a/src/cormaddp.f90 +++ /dev/null @@ -1,642 +0,0 @@ -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 quickselectdp, quickselectscalardp - end interface qselect - - interface iselect - procedure introselectdp, introselectscalardp - end interface iselect - - contains - - !!!!!!! DOUBLE PRECISION - 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 (atovector(subst) - end if - a => tovector(i+1) - if (atovector(subst) - end if - a => tovector(i+2) - if (atovector(subst) - end if - end do - do i=i,(r-1) - a => tovector(i) - if (atovector(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 (atovector(subst) - end if - a => tovector(i+1) - if (atovector(subst) - end if - a => tovector(i+2) - if (atovector(subst) - end if - end do - do i=i,(r-1) - a => tovector(i) - if (atovector(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 - !real(kind=dp) :: 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(int(i/5.0)+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(int(i/5.0)+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(int(i/5.0)+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(int(i/5.0)+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(int(i/5.0)+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 - - if (present(evencorrection) .or. present(ord)) then - continue - end if - - 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 - - if (present(evencorrection) .or. present(ord)) then - continue - end if - - quickselectscalardp=invector - end function quickselectscalardp - - -end module selectionalgo - -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 - !real(kind=dp) :: fresh - real(kind=dp), dimension(:), allocatable :: U, V - !real(kind=dp), dimension(:), allocatable :: 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 index 66209f7..32fd0f1 100644 --- a/src/init.c +++ b/src/init.c @@ -1,22 +1,66 @@ -#include -//#include +#include "RSCdefines.h" #include -#include -#include -#include #include #include +#include +#include +#include +#include +#include +#ifdef _OPENMP +#include +#endif + +// Wrapper for C cormad (with pointers) +SEXP cormad_C(SEXP R_matrix, SEXP len_rows, SEXP len_cols, SEXP R_evencorrect, + SEXP R_num_threads) { + // prepare input arguments + int n_row = asInteger(len_rows); + int n_col = asInteger(len_cols); + int output_size = (n_col - 1) * n_col / 2; + int evencorrect = + (n_row % 2 == 1) ? 0 : asInteger(R_evencorrect); // ingore for odd n + int num_threads = asInteger(R_num_threads); + SEXP output = PROTECT(allocVector(REALSXP, output_size)); // define output -extern void F77_NAME(cormadvecdp)(void *, void *, void *, void *, void *, void *); + // create pointers to (duplicate) matrix column + double *dupmat = + REAL(PROTECT(duplicate(R_matrix))); // copy input (to avoid modifications) + + // Call cormad +#ifdef _OPENMP + int max_threads = omp_get_max_threads(); + if (num_threads == 0) /* use max number of threads / 2 */ + num_threads = (max_threads / 2 == 0) ? 1 : max_threads / 2; + else if (num_threads > max_threads) + num_threads = max_threads; + else if (num_threads < 0) + num_threads = + (max_threads + num_threads > 0) ? max_threads + num_threads : 1; + + if (num_threads > 1) { + cormad_parallel(dupmat, n_row, n_col, REAL(output), evencorrect, + num_threads); + } else { + cormad(dupmat, n_row, n_col, REAL(output), evencorrect); + } +#else + if (num_threads != 1) + Rprintf("\nSpecified num_threads=%d, but OPENMP support not found; " + "switching to single core.\n\n", + num_threads); + cormad(dupmat, n_row, n_col, REAL(output), evencorrect); +#endif + // Remove protect to allow R grabage collect + UNPROTECT(2); + return output; +} -static const R_FortranMethodDef FortMethods[] = { - {"cormadvecdp", (DL_FUNC) &F77_NAME(cormadvecdp), 6}, - {NULL, NULL, 0} -}; +static const R_CallMethodDef callMethods[] = { + {"cormad_C", (DL_FUNC)&cormad_C, 5}, {NULL, NULL, 0}}; -void -R_init_RSC(DllInfo *dll) -{ - R_registerRoutines(dll, NULL, NULL, FortMethods, NULL); - R_useDynamicSymbols(dll, FALSE); +void R_init_RSC(DllInfo *dll) { + R_registerRoutines(dll, NULL, callMethods, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); + R_forceSymbols(dll, TRUE); } diff --git a/src/selection_algos_ptr.c b/src/selection_algos_ptr.c new file mode 100644 index 0000000..2abecac --- /dev/null +++ b/src/selection_algos_ptr.c @@ -0,0 +1,154 @@ +#include +/* Implements functions quickselect and introselect */ + +static double select_corner_cases(double *vector_shorter_3, int size_of_vector, + int extract_this_element_Ccorrected); + +static double pivot_of_3(double *vector_to_be_pivoted, int size); + +static int quickselect_onepass(double *vector, int size); + +/* Functions */ + +double quickselect_recursive(double *vector, int size, int k) { + /* Arguments: + * *vector: pointer to vector + * idN: size of vector + * k: k-th element to be extracted from the vector output: the + * extracted element of order k */ + + /* Decide if we need to go with quickselect_recursive or return output + * (quickselect_recursive works only for vectors with at least 3 elems) */ + if (size < 3) { // select the returning value + return select_corner_cases(vector, size, k); + } + + int subst = quickselect_onepass(vector, size); + + /* Determine where to go next (left of subst or right?) */ + if (subst == k) { + return *(vector + subst); + } else if (subst > k) { // go left + size = subst; + return quickselect_recursive(vector, size, k); + } else { // go right + // readjust k + k -= subst + 1; // take it back to orginal value (for full vector) + vector += subst + 1; + size -= subst + 1; + return quickselect_recursive(vector, size, k); + } +} + +static double select_corner_cases(double *vector, int size, int k) { + /* Used to handel corner cases not handeled by recursive strategy. + * Recursive strategy needs vectors of at least 3 elements + + * Arguments: + * *vector: pointer to vector of doubles to sort + * size: size of the vector + * k: k-th element to be extracted from the vector output: the + * extracted element of order k + */ + + // total elements in vector are idN-id0+1 + double ret = -111; + switch (size) { + case 1: + return *vector; + case 2: + switch (k) { + case 0: // return the smallest of the two + if (*vector < *(vector + 1)) + return *vector; + else + return *(vector + 1); + case 1: // return the biggest of the two + if (*vector > *(vector + 1)) + return *vector; + else + return *(vector + 1); + } + } + return ret; +} + +static double pivot_of_3(double *vector, int size) { + /* perform pivot of 3 returning the median element and + * rearranging the vector so to have: + * id0, i(see below) , idN -> min, max, med + + * Arguments + * vector: vector where to make substitutions + * id0, i, idN: positions of first, median and last element to consider + */ + + int idN = size - 1; // index of last element of the vector + int i = idN / 2; // index of middle element of the vector + double a, b, c; // auxiliary variables (avoid dereferncing too much) + a = *vector; + b = *(vector + i); + c = *(vector + idN); + + /* Nota; vogliamo che l'elemento mediano si trovi alla fine */ + double swapper; + if ((a > b) ^ (a > c)) { // id0 is median element + swapper = c; + c = a; + if (swapper > b) { + a = b; + b = swapper; + } else { + a = swapper; + } + } else if ((b > a) ^ (b > c)) { // i is median element + swapper = c; + c = b; + if (swapper > a) { + b = swapper; + } else { + b = a; + a = swapper; + } + } else { // idN is median element + if (a > b) { + swapper = a; + a = b; + b = swapper; + } + } + *vector = a; + *(vector + i) = b; + *(vector + idN) = c; + return c; +} + +static int quickselect_onepass(double *vector, int size) { + /* perform one pass for quickselect and returns the element where the pivot + * was substituted: this is used to decide whether to go left or right */ + + /* Pivoting section (pivot of 3) */ + double pivot = pivot_of_3(vector, size); + + /* One pass on vector */ + double swapper; + double *subst_el = vector; // address of element to substitute (copy not to + // alter main pointer) + int subst = 0; + for (int i = 0; i < size - 1; i++) { + if (*subst_el < pivot) { + if (subst == i) { + subst++; + } else { + swapper = *subst_el; + *subst_el = *(vector + subst); + *(vector + subst) = swapper; + subst++; + } + } + subst_el++; + } + *subst_el = *(vector + subst); + *(vector + subst) = pivot; + return subst; +} -- cgit v1.2.3