0.8.4: Cleaned comments and code

Author: slamballais

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: QDECR
 Type: Package
 Title: Vertex-wise statistical analysis
-Version: 0.8.3
+Version: 0.8.4
 Authors@R: as.person(c(
     "Ryan Muetzel <r.muetzel@erasmusmc.nl> [aut, cre]",
     "Sander Lamballais <s.lamballaistessensohn@erasmusmc.nl> [aut]"

NEWS.md

@@ -1,8 +1,8 @@
-# QDECR 0.8.3
+# QDECR 0.8.4
 
 ## Bug fixes
 * Fixed that the model output files (p.mgh, t.mgh, etc) also contained the mcz threshold (e.g. cache.th30).
-
+* Lots of comments removed/modified, and code was cleaned up a bit
 
 # QDECR 0.8.3
 

R/io_mgh.R

@@ -14,29 +14,25 @@ qdecr_prep_mgh <- function(input_path,
   if(measure2 == "w_g.pct") measure2 <- "w-g.pct"
   new_files <- file.path(input_path, files_list, "surf", paste(hemi, measure, fwhmc, target, "mgh", sep = "."))
   n <- length(new_files)
-
   temp <- load.mgh(new_files[1])
 
   m <- bigstatsr::FBM(temp$ndim1, n, backingfile = backing)
   cl <- if(verbose) parallel::makeForkCluster(n_cores, outfile = "") else parallel::makeForkCluster(n_cores)
   doParallel::registerDoParallel(cl)
-  
   capture.output(pb <- txtProgressBar(0, n, style = 3), file = "/dev/null")
-  
   foreach(i = seq_len(n)) %dopar% {
     setTxtProgressBar(pb, i)
     m[,i] <- load.mgh(new_files[i])$x
     NULL
   }
   parallel::stopCluster(cl)
-  
   if (length(mask) != nrow(m)) stop("Length of mask does not equal number of vertices from data.")
   m
 }
 
-
-
 #' Load an MGH file into memory
+#' 
+#' Originally written by Heath Perdoe (09/12/2013)
 #'
 #' @param input.file full path to the mgh file
 #'
@@ -45,7 +41,6 @@ qdecr_prep_mgh <- function(input_path,
 #' @export
 #' 
 load.mgh <- function(input.file) {
-  # written by Heath Pardoe, heath.pardoe at nyumc.org, 09/12/2013
   to.read <- file(input.file, "rb")
   v <- readBin(to.read, integer(), endian = "big")
   ndim1 <- readBin(to.read, integer(), endian = "big")
@@ -55,11 +50,8 @@ load.mgh <- function(input.file) {
   type <- readBin(to.read, integer(), endian = "big")
   dof <- readBin(to.read, integer(), endian = "big")
   close(to.read)
-
   to.read <- file(input.file,"rb")
   dump <- readBin(to.read, double(), size = 4, n = 71, endian = "big")
-  #THIS SEEMED TO ONLY READ IN THE FIRST CHUNK OF DATA, and ignored if additional subjects were merged in
-  #added ndim1*nframes
   x <- readBin(to.read ,double(), size = 4, n = ndim1*nframes, endian = "big")
   close(to.read)
   list(x = x, v = v, ndim1 = ndim1, ndim2 = ndim2, ndim3 = ndim3, nframes =
@@ -75,7 +67,6 @@ load.mgh <- function(input.file) {
 #' @export
 #' 
 load.annot <- function(input.file) {
-  # rewrite of load.mgh
   to.read <- file(input.file, "rb")
   on.exit (close(to.read))
   a_vtxct <- readBin(to.read, size = 4, integer(), endian = "big")
@@ -86,7 +77,6 @@ load.annot <- function(input.file) {
   a_len <- readBin(to.read, size = 4, integer(), endian = "big")
   a_fname <- readChar(to.read, a_len)
   a_num_entries <- readBin(to.read, size = 4, integer(), endian = "big")
-  
   LUT_Label <- LUT_len <- LUT_labelname <- LUT_red <- LUT_green <- LUT_blue <- LUT_transp <- rep(NA, a_num_entries)
   for (i in seq_len(a_num_entries)) {
     LUT_Label[i] <- readBin(to.read, size = 4, integer(), endian = "big")
@@ -99,11 +89,9 @@ load.annot <- function(input.file) {
   }
   close(to.read)
   on.exit()
-  
   a_vtxct2 <- seq(1, a_vtxct * 2, 2)
   a_vd_vno <- a_vd[a_vtxct2]
   a_vd_label <- a_vd[a_vtxct2 + 1]
-
   list(vtxct = a_vtxct, 
        vd_vno = a_vd_vno, 
        vd_label = a_vd_label,
@@ -127,18 +115,14 @@ load.annot <- function(input.file) {
 
 #' Save out an MGH file from a big matrix
 #'
+#' Modified version of save_mgh.m (by Heath Pardoe, 09/12/2013)
+#'
 #' @param vol MGH object (as from load.mgh)
 #' @param fname file name to be used to save out the data
 #'
 #' @export
 #' 
 bsfbm2mgh <-function(fbm, fname, filter = NULL) {
-
-  # R translation of save_mgh.m
-  # written by Heath Pardoe, heath.pardoe at nyumc.org, 09/12/2013
-  # modified by Ryan Muetzel to handle nframes/stacked data
-  # modified by Sander Lamballais to convert .bk files (from bigstatsr) to .mgh
-
   MRI.UCHAR <-  0
   MRI.INT <-    1
   MRI.LONG <-   2
@@ -156,24 +140,20 @@ bsfbm2mgh <-function(fbm, fname, filter = NULL) {
   } else {
     filter <- seq_len(fbm$nrow)
   }
-
   width <- fbm$ncol
   height <- 1
   depth <- 1
   nframes <- length(filter)
-
   writeBin(as.integer(1), fid, size = 4, endian = "big")
   writeBin(as.integer(width), fid, size = 4, endian = "big")
   writeBin(as.integer(height), fid, size = 4, endian = "big")
   writeBin(as.integer(depth), fid, size = 4, endian = "big")
   writeBin(as.integer(nframes), fid, size = 4, endian = "big")
   writeBin(as.integer(MRI.FLOAT), fid, size = 4, endian = "big")
   writeBin(as.integer(1), fid, size = 4, endian = "big")
-  
   UNUSED.SPACE.SIZE <- 256
   USED.SPACE.SIZE <- (3 * 4 + 4 * 3 * 4)
   unused.space.size <- UNUSED.SPACE.SIZE - 2
-
   writeBin(as.integer(0), fid, size = 2, endian = "big")
   writeBin(as.integer(rep.int(0, unused.space.size)), fid, size = 1)
   bpv <- 4
@@ -184,17 +164,15 @@ bsfbm2mgh <-function(fbm, fname, filter = NULL) {
 }
 
 #' Save out an MGH file from memory
+#' 
+#' #' R translation of save_mgh.m (by Heath Pardoe, 09/12/2013)
 #'
 #' @param vol MGH object (as from load.mgh)
 #' @param fname file name to be used to save out the data
 #'
 #' @export
 #'
 save.mgh <-function(vol,fname) {
-
-  # R translation of save_mgh.m
-  # written by Heath Pardoe, heath.pardoe at nyumc.org, 09/12/2013
-
   MRI.UCHAR <-  0
   MRI.INT <-    1
   MRI.LONG <-   2
@@ -203,90 +181,44 @@ save.mgh <-function(vol,fname) {
   MRI.BITMAP <- 5
   MRI.TENSOR <- 6
   slices <- c(1:256)
-
   fid <- file(fname, open = "wb", blocking = TRUE)
   on.exit(close(fid))
-
   width <- vol$ndim1
   height <- vol$ndim2
   depth <- vol$ndim3
-  #RLM added
   nframes <- vol$nframes
-
   writeBin(as.integer(1), fid, size = 4, endian = "big")
   writeBin(as.integer(width), fid, size = 4, endian = "big")
   writeBin(as.integer(height), fid, size = 4, endian = "big")
   writeBin(as.integer(depth), fid, size = 4, endian = "big")
-  #here we replace the default of 1 frame
   writeBin(as.integer(nframes), fid, size = 4, endian = "big")
-  #writeBin(as.integer(1),fid,size = 4, endian = "big")
-
-  # HP note: I've ignored all the diffusion tensor stuff
   writeBin(as.integer(MRI.FLOAT), fid, size = 4, endian = "big")
-
   writeBin(as.integer(1), fid, size = 4, endian = "big")
-  # dof = fread(fid, 1, 'int');
-  ## HP note: ignored ^this line^ from save_mgh.m
-
   UNUSED.SPACE.SIZE <- 256
-  USED.SPACE.SIZE <- (3 * 4 + 4 * 3 * 4)  # space for ras transform
-
+  USED.SPACE.SIZE <- (3 * 4 + 4 * 3 * 4)
   unused.space.size <- UNUSED.SPACE.SIZE - 2
-
-  # ignored all the stuff about "M" - could probably do it if necessary so let me know
-  #if (nargin > 2)
-  ## fwrite(fid, 1, 'short')        # ras.good.flag <- 0
-  # writeBin(1,fid,size = 2, endian = "big")
-  # unused.space.size <- unused.space.size - USED.SPACE.SIZE
-  ## fwrite(fid, sizes(1), 'float32')  # xsize
-  ## fwrite(fid, sizes(2), 'float32')  # ysize
-  ## fwrite(fid, sizes(3), 'float32')  # zsize
-  #
-  # fwrite(fid, M(1,1), 'float32')   # x.r
-  # fwrite(fid, M(2,1), 'float32')   # x.a
-  # fwrite(fid, M(3,1), 'float32')   # x.s
-
-  # fwrite(fid, M(1,2), 'float32')   # y.r
-  # fwrite(fid, M(2,2), 'float32')   # y.a
-  # fwrite(fid, M(3,2), 'float32')   # y.s
-
-  # fwrite(fid, M(1,3), 'float32')   # z.r
-  # fwrite(fid, M(2,3), 'float32')   # z.a
-  # fwrite(fid, M(3,3), 'float32')   # z.s
-
-  # fwrite(fid, M(1,4), 'float32')   # c.r
-  # fwrite(fid, M(2,4), 'float32')   # c.a
-  # fwrite(fid, M(3,4), 'float32')   # c.s
-  #else
-  # fwrite(fid, 0, 'short')        # ras.good.flag <- 0
   writeBin(as.integer(0), fid, size = 2, endian = "big")
-
-  #   } #
-
   writeBin(as.integer(rep.int(0, unused.space.size)), fid, size = 1)
-  bpv <- 4    # bytes/voxel
-  nelts <- width * height   # bytes per slice
-  #writeBin(vol$x, fid, size = 4, endian = "big")
+  bpv <- 4
+  nelts <- width * height
   writeBin(vol$x, fid, size = 4, endian = "big")
 }
 
-
 #' Create an object that is structured like an mgh object
 #'
 #' @param x the vertex-wise values
-#' @param v (to be added)
-#' @param ndim1 (to be added)
-#' @param ndim2 (to be added)
-#' @param ndim3 (to be added)
-#' @param nframes (to be added)
-#' @param type (to be added)
-#' @param dof (to be added)
+#' @param v Version (default = 1)
+#' @param ndim1 Width / 1st dimension
+#' @param ndim2 Height / 2nd dimension
+#' @param ndim3 Depth / 3rd dimension
+#' @param nframes Number of scalar components
+#' @param type Data type, can be UCHAR (0), SHORT (4), INT (1) or FLOAT (3)
+#' @param dof Degrees of freedom
 #'
 #' @export
 #'
 
 as_mgh <- function(x, v = NULL, ndim1 = NULL, ndim2 = NULL, ndim3 = NULL, nframes, type, dof) {
-  
   if (is.vector(x)) {
     v <- 1L
     ndim1 <- as.integer(length(x))
@@ -298,7 +230,6 @@ as_mgh <- function(x, v = NULL, ndim1 = NULL, ndim2 = NULL, ndim3 = NULL, nframe
   } else {
     stop("as_mgh only support objects of class `vector` right now.")
   }
-  
   out <- list(x = x, 
               v = v, 
               ndim1 = ndim1, 
@@ -309,5 +240,4 @@ as_mgh <- function(x, v = NULL, ndim1 = NULL, ndim2 = NULL, ndim3 = NULL, nframe
               dof = dof)
 
   return(out)
-  
 }

R/pool_quick.R

@@ -28,27 +28,21 @@ quick_pool <- function (qhat, se) {
   k <- length(qhat[[1]])
   im <- (1 + 1/m)
   um <- colMeans(do.call("rbind", se)^2)
-  
   qhat2 <- do.call("rbind", qhat)
   qbar <- colMeans(qhat2)
-  
   e <- sweep(qhat2, 2, qbar, `-`) 
   bm <- colSums(e^2)/(m - 1 + eps)
   t <- um + im * bm
   se2 <- sqrt(t)
   tval <- qbar/se2
-  
   lambda <- im * (bm/t)
-  
   eps2 <- 1e-04
   dfcom <- 1e+07
   lambda[lambda < eps2] <- eps2
   dfold <- (m - 1 + eps2)/lambda^2
   dfobs <- (dfcom + 1)/(dfcom + 3) * dfcom * (1 - lambda)
   df <- dfold * dfobs/(dfold + dfobs)
-  
   pval <- 2 * stats::pt(-abs(tval), df = df)
-  
   list(results = qbar, 
        se = se2, 
        t = tval, 
@@ -59,33 +53,26 @@ quick_pool2 <- function (qhat, se) {
 ### REWRITE OF mice::pool AND mice::mice_df to also handle multiple outcomes
   qhat <- lapply(qhat, as.matrix)
   se <- lapply(se, as.matrix)
-  
   eps <- 1e-100
   m <- length(qhat)
   dd <- dim(qhat[[1]])
   im <- (1 + 1/m)
   um <- Reduce("+", lapply(se, `^`, 2)) / m
-  
   qbar <- Reduce("+", qhat) / m 
   qhat2 <- array(unlist(qhat), dim = c(dd, m))
-  
   e <- sweep(qhat2, c(1,2),  qbar, `-`) 
   bm <- apply(e^2, c(1,2), sum)/(m - 1 + eps)
   t <- um + im * bm
   se2 <- sqrt(t)
   tval <- qbar/se2
-  
   lambda <- im * (bm/t)
-  
   eps2 <- 1e-04
   dfcom <- 1e+07
   lambda[lambda < eps2] <- eps2
   dfold <- (m - 1 + eps2)/lambda^2
   dfobs <- (dfcom + 1)/(dfcom + 3) * dfcom * (1 - lambda)
   df <- dfold * dfobs/(dfold + dfobs)
-  
   pval <- 2 * stats::pt(-abs(tval), df = df)
-  
   list(results = qbar, 
        se = se2, 
        t = tval,