R/bigcor_spiess.R

Defines functions anscore

Documented in anscore

#' A N Spiess bigcore function
#'
#' @export
anscore <- function( x, y = NULL, fun = c("cor", "cov"), size = 2000, verbose = TRUE, ...) {
    fun <- match.arg(fun)
  if (fun == "cor") FUN <- cor else FUN <- cov
    if (fun == "cor") STR <- "Correlation" else STR <- "Covariance" 
    if (!is.null(y) & NROW(x) != NROW(y)) stop("'x' and 'y' must have compatible dimensions!")
       
      NCOL <- ncol(x)
      if (!is.null(y)) YCOL <- NCOL(y)
          
        ## calculate remainder, largest 'size'-divisible integer and block size
        REST <- NCOL %% size
        LARGE <- NCOL - REST  
          NBLOCKS <- NCOL %/% size
            
          ## preallocate square matrix of dimension
          ## ncol(x) in 'ff' single format
          if (is.null(y)) resMAT <- ff(vmode = "double", dim = c(NCOL, NCOL))  
            else resMAT <- ff(vmode = "double", dim = c(NCOL, YCOL))
            
            ## split column numbers into 'nblocks' groups + remaining block
            GROUP <- rep(1:NBLOCKS, each = size)
              if (REST > 0) GROUP <- c(GROUP, rep(NBLOCKS + 1, REST))
            SPLIT <- split(1:NCOL, GROUP)
              
              ## create all unique combinations of blocks
              COMBS <- expand.grid(1:length(SPLIT), 1:length(SPLIT))
              COMBS <- t(apply(COMBS, 1, sort))
                COMBS <- unique(COMBS)  
                if (!is.null(y)) COMBS <- cbind(1:length(SPLIT), rep(1, length(SPLIT)))
                  
                  ## initiate time counter
                  timeINIT <- proc.time() 
                  
                  ## iterate through each block combination, calculate correlation matrix
                  ## between blocks and store them in the preallocated matrix on both
                  ## symmetric sides of the diagonal
                  for (i in 1:nrow(COMBS)) {
                        COMB <- COMBS[i, ]    
                      G1 <- SPLIT[[COMB[1]]]
                          G2 <- SPLIT[[COMB[2]]]    
                          
                          ## if y = NULL
                          if (is.null(y)) {
                                  if (verbose) cat(sprintf("#%d: %s of Block %s and Block %s (%s x %s) ... ", i, STR,  COMB[1],
                                                                                          COMB[2], length(G1),  length(G2)))      
                                RES <- FUN(x[, G1], x[, G2], ...)
                                      resMAT[G1, G2] <- RES
                                      resMAT[G2, G1] <- t(RES) 
                                          } else ## if y = smaller matrix or vector  
                                                {
                                                        if (verbose) cat(sprintf("#%d: %s of Block %s and 'y' (%s x %s) ... ", i, STR,  COMB[1],
                                                                                                                length(G1),  YCOL))    
                                            RES <- FUN(x[, G1], y, ...)
                                                  resMAT[G1, ] <- RES             
                                                }
                                          
                                          if (verbose) {
                                                  timeNOW <- proc.time() - timeINIT
                                                cat(timeNOW[3], "s\n")
                                                    }
                              
                              gc()
                                } 
                    
                    return(resMAT)
}
PietaSchofield/plibb documentation built on May 6, 2019, 6:45 p.m.