# R/corr2d.R In corr2D: Implementation of 2D Correlation Analysis in R

#### Documented in corr2dis.corr2d

#' Two-dimensional correlation analysis.
#'
#' \code{corr2d} calculates the synchronous and asynchronous correlation
#'     spectra between \code{Mat1} and \code{Mat1} (homo correlation)
#'     or between \code{Mat1} and \code{Mat2} (hetero correlation).
#'
#' \code{corr2d} uses a parallel fast Fourier transformation
#'     (\code{\link[stats]{fft}}) to calculate the complex correlation matrix.
#'     For parallelization the \code{\link[foreach]{foreach}} function is used.
#'     Large input matrices (> 4000 columns) can lead to long calculation times
#'     depending on the number of cores used. Also note that the resulting
#'     matrix can become very large, adjust the RAM limit with
#'     \code{\link[utils]{memory.limit}} accordingly. For a detailed description
#'     of the underlying math see references.
#'
#' @param Mat1,Mat2 Numeric matrix containing the data which will be correlated;
#'     '\emph{spectral variable}' by columns and '\emph{perturbation variables}'
#'     by rows. For hetero correlations \code{Mat1} and \code{Mat2} must have
#'     the same number of rows.
#' @param Ref1,Ref2 Numeric vector containing a single spectrum, which will be
#'     subtracted from \code{Mat1} (or \code{Mat2}, respectively) to generate dynamic spectra
#'     for 2D correlation analysis. Default is \code{NULL} in which case the \code{colMeans()}
#'     of \code{Mat1} (or \code{Mat2}, respectively) is used as reference. The length of \code{Ref1}
#'     (or \code{Ref2}) needs to be equal to the number of columns in \code{Mat1} (or \code{Mat2}).
#' @param Wave1,Wave2 Numeric vector containing the spectral variable. Needs to be
#'     specified if column names of \code{Mat1} (or \code{Mat2}) are undefined.
#' @param Time Numeric vector containing the perturbation variables. If specified, \code{Mat1}
#'     (and \code{Mat2} if given) will be interpolated to \code{N} equally spaced perturbation
#'     varibales using \code{Int} to speed up the fft algorithm.
#' @param Int Function specifying how the dataset will be interpolated to give
#'     \code{N} equally spaced perturbation variables. \code{\link[stats]{splinefun}}
#'     (default) or \code{\link[stats]{approxfun}} can for example be used.
#' @param N Positive, non-zero integer specifying how many equally spaced
#'     perturbation variables should be interpolated using \code{Int}. \code{N}
#'     should be higher than 4.\code{corr2d} is fastest if \code{N} is a power
#'     of 2.
#' @param Norm A number specifying how the correlation matrix should be
#'     normalized.
#' @param scaling Positive real number used as exponent when scaling the dataset
#'     with its standard deviation. Defaults to 0 meaning no scaling. 0.5
#'     (\emph{Pareto scaling}) and 1 (\emph{Pearson scaling}) are commonly used to enhance
#'     weak correlations relative to strong correlations.
#' @param corenumber Positive, non-zero integer specifying how many CPU cores
#'     should be used for parallel fft computation.
#' @param preview Logical: Should a 3D preview of the synchronous correlation
#'     spectrum be drawn at the end? Uses \code{\link[rgl]{persp3d}} from \pkg{rgl}
#'     package.
#'
#' @return \code{corr2D} returns a list of class "corr2d" containing the complex
#'     correlation matrix (\code{$FT}), the used reference spectra (\code{$Ref1},
#'     \code{$Ref2}), the spectral variables (\code{$Wave1}, \code{$Wave2}), the #' Fourier transformed data (\code{$ft1}, \code{$ft2}), the (interpolated) #' perturbation variables (\code{$Time}) and logical variable (\code{$Het}) #' indicating if homo (\code{FALSE}) or hetero (\code{TRUE}) correlation was done. #' #' @references #' I. Noda (1993) <DOI:10.1366/0003702934067694>\cr #' I. Noda (2012) <DOI:10.1016/j.vibspec.2012.01.006>\cr #' R. Geitner et al. (2019) <DOI:10.18637/jss.v090.i03> #' #' #' @seealso For plotting of the resulting list containing the 2D correlation #' spectra see \code{\link{plot_corr2d}} and \code{\link{plot_corr2din3d}}. #' #' @examples #' data(FuranMale, package = "corr2D") #' twod <- corr2d(FuranMale, Ref1 = FuranMale[1, ], corenumber = 1) #' #' plot_corr2d(twod, xlab = expression(paste("relative Wavenumber" / cm^-1)), #' ylab = expression(paste("relative Wavenumber" / cm^-1))) #' #' @aliases corr2d #' #' @export #' @importFrom foreach %dopar% #' @importFrom stats fft sd corr2d <- function(Mat1, Mat2 = NULL, Ref1 = NULL, Ref2 = NULL, Wave1 = NULL, Wave2 = NULL, Time = NULL, Int = stats::splinefun, N = 2^ceiling(log2(NROW(Mat1))), Norm = 1/(pi * (NROW(Mat1) - 1)), scaling = 0, corenumber = parallel::detectCores(), preview = FALSE) { # check user input for errors ----------------------------------------- if (!is.null(Ref1)) { if (!identical(length(Ref1), NCOL(Mat1))) { stop("length(Ref1) must be equal to NCOL(Mat1)") } } if (!is.null(Mat2) && !is.null(Ref2)) { if (!identical(length(Ref2), NCOL(Mat2))) { stop("length(Ref2) must be equal to NCOL(Mat2)") } } if (is.null(colnames(Mat1)) && is.null(Wave1)) { stop("Spectral variable 1 must be specified at Wave1 or in colnames(Mat1)") } if (!is.null(Mat2)) { if (is.null(colnames(Mat2)) && is.null(Wave2)) { stop("Spectral variable 2 must be specified at Wave2 or in colnames(Mat2)") } } if (N <= 0 || N%%1 != 0) { stop("N must be a positive, non-zero integer") } if (corenumber <= 0 || corenumber%%1 != 0) { stop("corenumber must be a positive, non-zero integer") } if (!is.logical(preview)) { stop("preview needs to be logical") } # create an R session for every detected CPU core for parallel computing cl <- parallel::makeCluster(corenumber) doParallel::registerDoParallel(cl) # define Wave1/2 from colnames if needed ------------------------------ if (is.null(Wave1)) { Wave1 <- as.numeric(colnames(Mat1)) } if (!is.null(Mat2) && !identical(Mat1, Mat2) && is.null(Wave2)) { Wave2 <- as.numeric(colnames(Mat2)) } # spectral variable must be increasing, also switch Ref if needed ------ if (is.unsorted(Wave1)) { Wave1 <- Wave1[length(Wave1):1] Mat1 <- Mat1[, NCOL(Mat1):1] if (!is.null(Ref1)) { Ref1 <- Ref1[length(Ref1):1] } } if (!is.null(Mat2) && !identical(Mat1, Mat2) && is.unsorted(Wave2)) { Wave2 <- Wave2[length(Wave2):1] Mat2 <- Mat2[, NCOL(Mat2):1] if (!is.null(Ref2)) { Ref2 <- Ref2[length(Ref2):1] } } # define Wave2 (if needed), D1/2, Het --------------------------------- if (is.null(Mat2) || identical(Mat1, Mat2)) { cat(c("HOMO-Correlation:", corenumber, "cores used for calculation\n")) Het <- FALSE D1 <- D2 <- dim(Mat1) Wave2 <- Wave1 Mat2 <- Mat1 Ref2 <- Ref1 } else { cat(c("HETERO-Correlation:", corenumber, "cores used for calculation\n")) Het <- TRUE D1 <- dim(Mat1) D2 <- dim(Mat2) } # Interpolate values of perturbation variable for equidistant --------- # datapoints for fft -------------------------------------------------- if ((is.null(Time) == FALSE) & (length(Time) == D1[1])) { cat(c(format(Sys.time(), "%X -"), "Interpolate Time from", min(Time), "to", max(Time), "\n", "to obtain", N, "equidistant data points for FFT", "\n")) TIME <- seq(min(Time), max(Time), length.out = N) tmp <- apply(Mat1, 2, function(y) Int(x = Time, y = y)) Mat1 <- sapply(tmp, function(x) x(TIME)) if (Het == FALSE) { Mat2 <- Mat1 } else { tmp <- apply(Mat2, 2, function(y) Int(x = Time, y = y)) Mat2 <- sapply(tmp, function(x) x(TIME)) } Time <- TIME } # Get perturbation variables from rownames ---------------------------- if (is.null(Time) && !is.null(rownames(Mat1))) { Time <- as.numeric(rownames(Mat1)) } # Subtract reference ------------------------------------------------- if (is.null(Ref1)) { cat(c(format(Sys.time(), "%X -"), "using mean values as reference\n")) Ref1 <- colMeans(Mat1) } Mat1 <- sweep(Mat1, 2, Ref1, "-") if (Het == FALSE) { Mat2 <- Mat1 Ref2 <- Ref1 } else { if (is.null(Ref2)) { Ref2 <- colMeans(Mat2) } Mat2 <- sweep(Mat2, 2, Ref2, "-") } # Apply scaling ------------------------------------------------------- if (scaling > 0) { cat(c(format(Sys.time(), "%X -"), "apply scaling\n")) sd1 <- apply(Mat1, 1, sd) sd2 <- apply(Mat2, 1, sd) Mat1 <- Mat1 / (sd1^scaling) Mat2 <- Mat2 / (sd2^scaling) } # Do fft for every wavenumber in both matrices to obtain frequency ---- # dependent vectors for every wavenumber and scalar-multiply them ----- # for every point in the 2D-correlation spectrum ---------------------- i <- NULL # Dummy to trick R CMD check FT <- NULL cat(c(format(Sys.time(), "%X -"), "Fast Fourier Transformation and multiplication \n", "to obtain a", D1[2], "x", D2[2], "correlation matrix", "\n")) FT <- matrix(NA, NCOL(Mat1), NCOL(Mat2)) # Selection of fft elements ftele1 <- if(NROW(Mat1) %% 2 == 0) { 1:(NROW(Mat1) - 1) %/% 2 + 1 } else { 1:(NROW(Mat1)) %/% 2 + 1 } ft1 <- foreach::foreach(i = 1:NCOL(Mat1), .combine = 'cbind') %dopar% { fft(Mat1[, i])[ftele1] } if (Het == FALSE) { ft2 <- ft1 } else { # Selection of fft elements ftele2 <- if(NROW(Mat2) %% 2 == 0) { 1:(NROW(Mat2) - 1) %/% 2 + 1 } else { 1:(NROW(Mat2)) %/% 2 + 1 } ft2 <- foreach::foreach(i = 1:NCOL(Mat2), .combine = 'cbind') %dopar% { fft(Mat2[, i])[ftele2] } } FT <- matrix(Norm * parallel::parCapply(cl, ft1, get("%*%"), Conj(ft2)), NCOL(ft1), NCOL(ft2), byrow = TRUE) cat(c(format(Sys.time(), "%X -"), "Done\n")) Obj <- list(FT = FT, Ref1 = Ref1, Ref2 = Ref2, Wave1 = Wave1, Wave2 = Wave2, ft1 = ft1, ft2 = ft2, Time = Time, Het = Het) parallel::stopCluster(cl) # 3d preview of the synchronous correlation spectrum ------------------ if (preview == TRUE) { if (dim(FT)[1] > 700) { tmp1 <- round(seq(1, dim(FT)[1], length = 700)) } else { tmp1 <- seq(1, dim(FT)[1], 1) } if (dim(FT)[2] > 700) { tmp2 <- round(seq(1, dim(FT)[2], length = 700)) } else { tmp2 <- seq(1, dim(FT)[2], 1) } rgl::persp3d(Wave1[tmp1], Wave2[tmp2], Re(FT)[tmp1, tmp2], col = "grey") } class(Obj) <- "corr2d" return(Obj) } #' @method summary corr2d #' @export summary.corr2d <- function(object, ...) { cat(c(NROW(object$FT), "x", NCOL(object$FT), if (object$Het == TRUE){"Hetero"} else {"Homo"}, "correlation spectra\n"))
if (object$Het == TRUE) { cat(c("Spectral variable 1:", min(object$Wave1), "-", max(object$Wave1), "\n")) cat(c("Spectral variable 2:", min(object$Wave2), "-", max(object$Wave2), "\n")) } else { cat(c("Spectral variable:", min(object$Wave1), "-", max(object$Wave1), "\n")) } if (!is.null(object$Time)) {
cat(c("Perturbation variable:", length(object$Time), "values,", min(object$Time), "-", max(object\$Time), "\n"))
}

}

#' @method print corr2d
#' @export
print.corr2d <- function(x, ...)
{
print.default(x, ...)

}

#' Check for object class "corr2d"
#'
#' The function checks if an object is of class "corr2d".
#'
#' The function uses the \code{\link{inherits}} function.
#'
#' @param x An object which should be check if it is of class "corr2d".
#'
#' @return A logical scalar
#'
#' @examples
#'     data(FuranMale, package = "corr2D")
#'     twod <- corr2d(FuranMale, Ref1 = FuranMale[1, ], corenumber = 1)
#'
#'     # TRUE
#'     is.corr2d(twod)
#'     # FALSE
#'     is.corr2d(2)
#'
#' @references
#'     R. Geitner et al. (2019) <DOI:10.18637/jss.v090.i03>
#'
#' @export
is.corr2d <- function(x)
{
inherits(x, "corr2d", which = FALSE)
}


## Try the corr2D package in your browser

Any scripts or data that you put into this service are public.

corr2D documentation built on July 14, 2022, 5:06 p.m.