reimplemented prep_fun; something had gone wrong when reverting a branch

Author: slamballais

R/qdecr.R

@@ -104,12 +104,12 @@ qdecr <- function(id,
   md <- imp2list(data)
 
   # Assemble and check input arguments
-  vw$input <- qdecr_check(id, md, margs, hemi, vertex, measure, model, target, project, dir_out_tree, clobber, fwhm, n_cores)
+  vw$input <- qdecr_check(id, md, margs, hemi, vertex, measure, model, target, project, dir_out_tree, clobber, fwhm, n_cores, prep_fun)
   vw$mask <- qdecr_check_mask(mask, mask_path)
   vw$paths <- check_paths(vw, dir_tmp, dir_subj, dir_out, dir_fshome, mask_path)
 
   # Check model
-  vw$model <- qdecr_model(vw$input$model, vw$input$md, vw$input$id, vw$input$vertex, vw$input$margs, vw$paths$dir_tmp2)
+  vw$model <- qdecr_model(vw$input$model, vw$input$prep_fun, vw$input$md, vw$input$id, vw$input$vertex, vw$input$margs, vw$paths$dir_tmp2)
 
   # Check backing
   qdecr_check_backing(c(vw$model$backing, vw$model$backing_to_remove, vw$paths$backing_mgh), clobber)

R/qdecr_check.R

@@ -14,7 +14,7 @@ qdecr_check_backing <- function(backing, clobber){
 }
 
 qdecr_check <- function(id, md, margs, hemi, vertex, measure, model, target,
-                        project, dir_out_tree, clobber, fwhm, n_cores){
+                        project, dir_out_tree, clobber, fwhm, n_cores, prep_fun){
 
   # verify that id is correct and exists within data
   check_id(id, md)
@@ -27,6 +27,9 @@ qdecr_check <- function(id, md, margs, hemi, vertex, measure, model, target,
 
   # Check whether the number of cores is correct
   check_cores(n_cores)
+  
+  # Check whether prep_fun exists
+  check_prep_fun(prep_fun)
 
   # Return all output arguments assembled
   input <- list()
@@ -45,6 +48,7 @@ qdecr_check <- function(id, md, margs, hemi, vertex, measure, model, target,
   input[["fwhm"]] <- fwhm
   input[["fwhmc"]] <- paste0("fwhm", fwhm)
   input[["n_cores"]] <- n_cores
+  input[["prep_fun"]] <- prep_fun
   input
 }
 
@@ -167,6 +171,11 @@ check_paths <- function(vw, dir_tmp, dir_subj, dir_out, dir_fshome, mask_path){
  paths
 }
 
+check_prep_fun <- function(prep_fun){
+  tryCatch(get2(prep_fun), error = function(e) stop("Provided `prep_fun` cannot be found."))
+  invisible(NULL)
+}
+
 check_vertex <- function(vertex, id, md){
   stxt <- c("Please provide a unique, unused name to `vertex`. \n",
             "We only use it as an identifier in the model output.")

R/qdecr_model.R

@@ -1,44 +1,22 @@
-
-qdecr_model <- function(model, md, id, vertex, margs, dir_tmp2){
-  
+qdecr_model <- function(model, prep_fun, md, id, vertex, margs, dir_tmp2) {
   prepvw <- list(model = model, data = md, id = id, 
-             vertex = vertex, margs = margs)
+                 vertex = vertex, margs = margs, path = dir_tmp2)
   class(prepvw) <- "prepvw"
-  
-  if (prepvw$model == "QDECR::megha") {
-    model_qdecr_megha(prepvw)
-  } else if (prepvw$model == "RcppEigen::fastLm") {
-      prepvw$so <- c("coef", "se", "t", "p", "resid")
-      prepvw$backing <- paste0(dir_tmp2, "_", prepvw$so[1:4], "_backend")
-      prepvw$backing_to_remove <- paste0(dir_tmp2, "_", prepvw$so[5], "_backend")
-      model_RcppEigen_fastLm(prepvw) 
-  } else if (prepvw$model == "stats::glm") {
-      model_stats_glm(prepvw)
-  } else if (prepvw$model == "stats::lm") {
-      model_stats_lm(prepvw)
-  } else if (prepvw$model == "survival::coxph") {
-      model_survival_coxph(prepvw)
-  } else if (prepvw$model == "default") {
-      prepvw$so <- c("coef", "se", "t", "p", "resid")
-      prepvw$backing <- paste0(dir_tmp2, "_", prepvw$so[1:4], "_backend")
-      prepvw$backing_to_remove <- paste0(dir_tmp2, "_", prepvw$so[5], "_backend")
-      model_default(prepvw)
-  }
+  do.call2(prep_fun, list(prepvw = prepvw))
 }
 
-model_default <- function(prepvw) {
-  prepvw$ff <- "vw_default"
-  prepvw$formula <- prepvw$margs$formula
-  vw <- prepvw
-  class(vw) <- c(vw$ff, "vw")
-  vw
+qdecr_backing_path <- function(prepvw, to_keep, to_remove) {
+  prepvw$so <- c(to_keep, to_remove)
+  prepvw$backing <- paste0(prepvw$path, "_", prepvw$so[1:4], "_backend")
+  prepvw$backing_to_remove <- paste0(prepvw$path, "_", prepvw$so[5], "_backend")
+  return(prepvw)
 }
 
-model_qdecr_megha <- function(prepvw){
-  stop("`QDECR::megha` has not been implemented yet.")
-}
-
-model_RcppEigen_fastLm <- function(prepvw){
+prep_fastlm <- function(prepvw){
+  to_keep <- c("coef", "se", "t", "p")
+  to_remove <- "resid"
+  prepvw <- qdecr_backing_path(prepvw, to_keep, to_remove)
+  
   mf <- prepvw$margs
   if (is.null(mf$formula)) 
     stop("No `formula` set in margs.")
@@ -63,15 +41,15 @@ model_RcppEigen_fastLm <- function(prepvw){
   mx_test[, prepvw$vertex] <- 999
 
   if (!identical(mx[[nn]], mx_test)) stop ("Somewhere in your formula you specified a special term related to your vertex measure", 
-    " (interaction, polynomial, AsIs, etc); `qdecr_fastlm` currently does not support this.")
-
+                                           " (interaction, polynomial, AsIs, etc); `qdecr_fastlm` currently does not support this.")
+  
   y <- model.response(mx[[nn]], "numeric")
 
   if (nrow(mx[[nn]]) != nr) stop("The data that you are putting into the regression has missings! \n",
                                  "QDECR can't handle that yet; we will fix this soon!")
 
   ys <- if(identical(unname(y), rep(999, nrow(mx[[nn]])))) "LHS" else "RHS"
-   
+  
   if (ys == "LHS") { 
     mx_test2 <- mx_test
     mx_test2b <- model.matrix(mx_test2, object = mt, contrasts)
@@ -92,78 +70,5 @@ model_RcppEigen_fastLm <- function(prepvw){
   }
   class(vw) <- c(vw$ff, "vw")
   vw
+  
 }
-
-model_stats_glm <- function(prepvw){
-  stop("`stats::glm` has not been implemented yet.")
-}
-
-model_stats_lm <- function(prepvw){
-  mf <- prepvw$margs
-  ret.x <- mf$x
-  ret.y <- mf$y
-  if (is.null(mf$drop.unused.levels)) 
-    mf$drop.unused.levels <- TRUE
-  if (is.null(mf$formula)) 
-    stop("No `formula` set in margs.")
-  iii <- c("formula", "data", "subset", "weights", "na.action", "offset")
-  mf[iii] <- mf[!sapply(mf, is.symbol)][iii]
-  mfz <- mf[match(c("formula", "data", "subset", "weights", "na.action", "offset"), names(mf), 0L)]
-  nr <- nrow(prepvw$data[[1]])
-  mx <- lapply(prepvw$data, function(x) {
-    mfz$data <- x
-    mfz$data[, prepvw$vertex] <- 999
-    do.call2("stats::model.frame", mfz)
-  })
-  nn <- length(prepvw$data)
-  if (mf$method == "model.frame")
-    stop("`margs$method` is equal to `model.frame` (for the lm call). Try `qr`.") else 
-      if (mf$method != "qr")
-    stop(gettextf("margs$method = '%s' is not supported in lm. Try 'qr'",
-                     method), domain = NA)
-  mt <- attr(mx[[nn]], "terms")
-  y <- model.response(mx[[nn]], "numeric")
-  ys <- if(identical(unname(y), rep(999, nr))) "LHS" else "RHS"
-  w <- as.vector(model.weights(mx[[nn]]))
-  if (!is.null(w) && !is.numeric(w)) 
-    stop("'margs$weights' must be a numeric vector (for the lm call)")
-  offset <- as.vector(model.offset(mf))
-  if (!is.null(offset)) {
-    if (length(offset) != NROW(y)) 
-      stop(gettextf("number of offsets is %d, should equal %d (number of observations) (in lm call)", 
-                    length(offset), NROW(y)), domain = NA)
-  }
-  if (is.empty.model(mt)) stop("The provided model (to lm) is empty. Check your data + formula.")
-  test <- model.matrix(mt, mx[[1]], contrasts) 
-  test2 <- model.matrix(mt, mx[[nn]], contrasts)
-  mm <- NULL
-  prepvw$ff <- "vw_lm_fit_slow"
-  if (prepvw$vertex %in% colnames(attr(mt, "factors")) || ys == "LHS"){
-    # mfz2 <- mfz
-    # mfz2$data <- prepvw$data[[1]]
-    # mfz2$data[, prepvw$vertex] <- mfz$data[, prepvw$vertex]
-    # mm2 <- do.call2("stats::model.frame", mfz2)
-    if(identical(test, test2)){
-      mm <- lapply(mx, model.matrix, object = mt, contrasts)
-      ff <- if (is.null(w)) "vw_lm_fit" else "vw.lm_wfit"
-      vw <- list(mm = mm, ff = ff, vertex = prepvw$vertex, ys = ys,
-                 y = y, w = w, offset = offset,
-                 singular.ok = prepvw$margs$singular.ok)
-    } else {
-      warning("Your formula for `lm` contains computed terms. \n", 
-              "We will rely on the slower implementation of our lm.")
-      vw <- prepvw
-    }
-  } else {
-    warning("Your formula for `lm` contains complicated terms. \n",
-            "We will rely on the slower implementation of our lm.")
-    vw <- prepvw
-  }
-  class(vw) <- c(vw$ff, "vw")
-  vw
-}
-
-model_survival_coxph <- function(prepvw){
-  stop("`survival::coxph` has not been implemented yet.")
-}
-